Install Packages

#library(nlme)
#library(ggplot2)
#library(dplyr)
#library(broom)
#library(purrr)
#library(tidyr)
#library(readxl)
#library(RColorBrewer)
#library(gridExtra)
#library(vegan)
#library(patchwork)
#library(cowplot)
#library(grid)
#library(lme4)
#library(plotly)
#library(MuMIn)
#library(lmerTest)
#library(arm)
#library(sjPlot)
#library(indicspecies)
#library(stringr)

Load Packages

library(nlme)
library(ggplot2)
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following object is masked from 'package:nlme':
## 
##     collapse
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(broom)
library(purrr)
library(tidyr)
library(readxl)
library(RColorBrewer)
library(gridExtra)
## 
## Attaching package: 'gridExtra'
## The following object is masked from 'package:dplyr':
## 
##     combine
library(vegan)
## Loading required package: permute
## Loading required package: lattice
library(patchwork)
library(cowplot)
## 
## Attaching package: 'cowplot'
## The following object is masked from 'package:patchwork':
## 
##     align_plots
library(grid)
library(lme4)
## Loading required package: Matrix
## 
## Attaching package: 'Matrix'
## The following objects are masked from 'package:tidyr':
## 
##     expand, pack, unpack
## 
## Attaching package: 'lme4'
## The following object is masked from 'package:nlme':
## 
##     lmList
library(plotly)
## 
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## The following object is masked from 'package:stats':
## 
##     filter
## The following object is masked from 'package:graphics':
## 
##     layout
library(MuMIn)
## Registered S3 methods overwritten by 'MuMIn':
##   method        from 
##   nobs.multinom broom
##   nobs.fitdistr broom
library(lmerTest)
## 
## Attaching package: 'lmerTest'
## The following object is masked from 'package:lme4':
## 
##     lmer
## The following object is masked from 'package:stats':
## 
##     step
library(arm)
## Loading required package: MASS
## 
## Attaching package: 'MASS'
## The following object is masked from 'package:plotly':
## 
##     select
## The following object is masked from 'package:patchwork':
## 
##     area
## The following object is masked from 'package:dplyr':
## 
##     select
## 
## arm (Version 1.14-4, built: 2024-4-1)
## Working directory is /cloud/project
## 
## Attaching package: 'arm'
## The following object is masked from 'package:MuMIn':
## 
##     coefplot
library(sjPlot)
## Install package "strengejacke" from GitHub (`devtools::install_github("strengejacke/strengejacke")`) to load all sj-packages at once!
## 
## Attaching package: 'sjPlot'
## The following objects are masked from 'package:cowplot':
## 
##     plot_grid, save_plot
library(indicspecies)
library(stringr)

Functions

save_plot <- function(plot_obj, file_name, width = 8, height = 6, dpi = 300, format = "png", plot_title = NULL) {
  format <- tolower(format)

  # Update plot title if provided
  if (!is.null(plot_title)) {
    plot_obj <- plot_obj + ggtitle(plot_title)
  }

  # Choose graphics device
  if (format == "png") {
    png(paste0(file_name, ".png"), width = width, height = height, units = "in", res = dpi)
  } else if (format == "pdf") {
    pdf(paste0(file_name, ".pdf"), width = width, height = height)
  } else if (format == "jpeg") {
    jpeg(paste0(file_name, ".jpeg"), width = width, height = height, units = "in", res = dpi)
  } else {
    stop("Unsupported file format")
  }

  # Print and close
  print(plot_obj)
  dev.off()
}



make_seedling_plot <- function(data, group_label, color_line = "steelblue", color_point = "steelblue4", custom_title = NULL) {
  
  # Fit model
  model <- lme(fixed = sum_seedlings ~ Carbon_tonne_ha, 
               random = ~ 1 | Site, 
               data = data)
  
  # Summary stats
  model_summary <- summary(model)
  p_val <- model_summary$tTable["Carbon_tonne_ha", "p-value"]
  slope <- model_summary$tTable["Carbon_tonne_ha", "Value"]
  
  r2_vals <- r.squaredGLMM(model)
  r2_marginal <- r2_vals[1]
  r2_conditional <- r2_vals[2]
  
  label_text <- paste0(
    "beta = ", round(slope, 3), "\n",
    "R2[m] = ", round(r2_marginal, 3), "\n",
    "R2[c] = ", round(r2_conditional, 3), "\n",
    "p = ", signif(p_val, 3)
  )
  
  # Plot
  graph <- plot_model(model, type = "pred", colors = color_line) +
    geom_point(data = data, aes(x = Carbon_tonne_ha, y = sum_seedlings), 
               color = color_point, size = 1.5) +
    labs(
      x = "Carbon (tonne/ha)",
      y = "Total Seedlings",
      title = ifelse(is.null(custom_title), paste("Seedling Abundance vs topsoil SOC (", group_label, ")", sep = ""), custom_title)
    ) +
    annotate("text", x = Inf, y = Inf, label = label_text, 
             hjust = 1.1, vjust = 1.5, size = 4, family = "mono") +
    theme_bw() +
    theme(
      )
    
  
  return(graph)
}

Import Data

carbon_dataset <- read_excel("/cloud/project/Remnant soil communities - data - Mallie - Loyne - Cluanie + inverts + carbon.xlsx", sheet = 10)

overview_dataset <- read_excel("/cloud/project/Remnant soil communities - data - Mallie - Loyne - Cluanie + inverts + carbon.xlsx", sheet = 2)


seedling_dataset <- read_excel("/cloud/project/Remnant soil communities - data - Mallie - Loyne - Cluanie + inverts + carbon.xlsx", sheet = 6)

depth_dataset <- read_excel("/cloud/project/Remnant soil communities - data - Mallie - Loyne - Cluanie + inverts + carbon.xlsx", sheet = 9)
## New names:
## • `` -> `...10`
plant_com_dataset <- read_excel("/cloud/project/Remnant soil communities - data - Mallie - Loyne - Cluanie + inverts + carbon.xlsx", sheet = 3)

pitfall_dataset <- read_excel("/cloud/project/Remnant soil communities - data - Mallie - Loyne - Cluanie + inverts + carbon.xlsx", sheet = 8)

soil_invert_dataset <- read_excel("/cloud/project/Remnant soil communities - data - Mallie - Loyne - Cluanie + inverts + carbon.xlsx", sheet = 11)

Process Data

#Mutate overview dataset to convert "Unforested (heath)",  "Unforested (heathland)", "Unforested (heath)" and "Unforested (remnant_scat)" to "Unforested (remnant)"

overview_dataset <- overview_dataset %>%
  mutate(Type = ifelse(Type %in% c("Unforested (heath)", "Unforested (heathland)"), "Unforested (heath)", Type)) %>%
  mutate(Type = ifelse(Type %in% c("Unforested (remnant_scat)"), "Unforested (remnant)", Type))

carbon_dataset <- merge(carbon_dataset, overview_dataset, by = "Plot")

plant_com_dataset <- merge(plant_com_dataset, overview_dataset, by = "Plot")


plant_com_dataset <- plant_com_dataset %>%
  mutate(Species = str_replace(Species, "Clubmoss", "Club moss"))

# Subset the data set to include only relevant information:

carbon_dataset <- carbon_dataset[, c("Site", "Plot","Layer" ,"Type", "Carbon (tonne/ha)")]

carbon_dataset <- carbon_dataset %>%
  rename(Carbon_tonne_ha = `Carbon (tonne/ha)`)


# Subset the data set to include only relevant information:

carbon_dataset <- carbon_dataset[, c("Site", "Plot","Layer" ,"Type", "Carbon_tonne_ha")]


sum_seedlings <- seedling_dataset %>%
  group_by(Plot) %>%
  summarise(sum_seedlings = sum(Frequency), na.rm = TRUE)

est_seedlings <- seedling_dataset %>%
  filter(Height_category %in% c("100-150", ">150", "50-100"))

est_seedlings <- est_seedlings %>%
  group_by(Plot) %>%             # Group by Plot
  summarise(Total_Seedlings = sum(Frequency, na.rm = TRUE))

est_seedlings$sum_seedlings <- NULL                      # Remove the column
names(est_seedlings)[names(est_seedlings) == "Total_Seedlings"] <- "sum_seedlings"  # Rename

carbon_dataset <- carbon_dataset %>%
  filter(`Carbon_tonne_ha` != 0)

# First, get the plots with 0-10, O
o_samples <- carbon_dataset %>%
  filter(Layer == "0-10, O")

# Then, get the plots with 0-10, M, but only those that do not have 0-10, O
m_samples <- carbon_dataset %>%
  filter(Layer == "0-10, M") %>%
  anti_join(o_samples, by = "Plot")  # Exclude plots that have 0-10, O

# Combine both datasets (0-10, O + the filtered 0-10, M)
upper10_data <- bind_rows(o_samples, m_samples)

upper10_seedlings <- merge(upper10_data, sum_seedlings, by = "Plot")

upper10_est_seedlings <- merge(upper10_data, est_seedlings, by = "Plot")

regen_upper10 <- upper10_seedlings[upper10_seedlings$Type == "Regen", ]

regen_upper10_est <- upper10_est_seedlings[upper10_est_seedlings$Type == "Regen", ]

carbon_0_20 <- carbon_dataset %>% 
  # Create a new layer that combines "0-10, O" and "10-20, O"
  mutate(Layer_combined = case_when(
    Layer == "0-10, O" ~ "0-20, O",
    Layer == "10-20, O" ~ "0-20, O",
    TRUE ~ Layer # Keep other layers unchanged
  )) %>%
  # Now group by the combined layer and other factors (Type, Site)
  group_by(Plot, Type, Layer_combined, Site) %>%
  summarise(
    Total_Carbon_Stock = sum(Carbon_tonne_ha) # Upper bound of 95% CI
  )
## `summarise()` has grouped output by 'Plot', 'Type', 'Layer_combined'. You can
## override using the `.groups` argument.
plant_com_dataset <- plant_com_dataset %>%
  mutate(DOMIN = case_when(
    DOMIN == 1 ~ 1,
    DOMIN == 2 ~ 2,
    DOMIN == 3 ~ 3,
    DOMIN == 4 ~ 7,
    DOMIN == 5 ~ 18,
    DOMIN == 6 ~ 29.5,
    DOMIN == 7 ~ 42,
    DOMIN == 8 ~ 63,
    DOMIN == 9 ~ 83,
    DOMIN == 10 ~ 95.5  # For any values outside 1-10
  ))

plant_com_dataset <- plant_com_dataset[c("Plot", "Site.y", "Species", "DOMIN", "Type")]

plant_rich <- plant_com_dataset %>%
  group_by(Plot) %>%
  summarise(Richness = sum(DOMIN > 0))

plant_rich <- merge(plant_rich, overview_dataset, by = "Plot")

# Calculate Simpson's Diversity Index for each plot
plant_div <- plant_com_dataset %>%
  group_by(Plot) %>%
  # Calculate total abundance for each plot
  mutate(Total_DOMIN = sum(DOMIN)) %>%
  # Calculate relative abundance for each species
  mutate(Relative_DOMIN = DOMIN / Total_DOMIN) %>%
  # Calculate Simpson's D (sum of squared relative abundances)
  summarise(Simpsons_D = sum(Relative_DOMIN^2, na.rm = TRUE)) %>%
  # Calculate Simpson's Diversity Index (1 - D)
  mutate(Diversity_Index = 1 - Simpsons_D)

plant_div <- merge(plant_div, overview_dataset, by = "Plot")


soil_invert_dataset <- merge(soil_invert_dataset, overview_dataset, by = "Plot")

soil_invert_dataset <- soil_invert_dataset[c("Plot", "Site.y", "Order", "Abundance", "Type")]

soil_invert_rich <- soil_invert_dataset %>%
  group_by(Plot) %>%
  summarise(Richness = sum(Abundance > 0))

soil_invert_rich <- merge(soil_invert_rich, overview_dataset, by = "Plot")

# Calculate Simpson's Diversity Index for each plot
soil_invert_div <- soil_invert_dataset %>%
  group_by(Plot) %>%
  # Calculate total abundance for each plot
  mutate(Total_Abundance = sum(Abundance)) %>%
  # Calculate relative abundance for each species
  mutate(Relative_Abundance = Abundance / Total_Abundance) %>%
  # Calculate Simpson's D (sum of squared relative abundances)
  summarise(Simpsons_D = sum(Relative_Abundance^2, na.rm = TRUE)) %>%
  # Calculate Simpson's Diversity Index (1 - D)
  mutate(Diversity_Index = 1 - Simpsons_D)

soil_invert_div <- merge(soil_invert_div, overview_dataset, by = "Plot")

pitfall_dataset <- merge(pitfall_dataset, overview_dataset, by = "Plot")

pitfall_dataset <- pitfall_dataset[c("Plot", "Site.y", "Order", "Abundance", "Type")]

pitfall_rich <- pitfall_dataset %>%
  group_by(Plot) %>%
  summarise(Richness = sum(Abundance > 0))

pitfall_rich <- pitfall_rich %>%
  filter(Richness > 0)

pitfall_rich <- merge(pitfall_rich, overview_dataset, by = "Plot")

# Calculate Simpson's Diversity Index for each plot
pitfall_div <- pitfall_dataset %>%
  group_by(Plot) %>%
  # Calculate total abundance for each plot
  mutate(Total_Abundance = sum(Abundance)) %>%
  # Calculate relative abundance for each species
  mutate(Relative_Abundance = Abundance / Total_Abundance) %>%
  # Calculate Simpson's D (sum of squared relative abundances)
  summarise(Simpsons_D = sum(Relative_Abundance^2, na.rm = TRUE)) %>%
  # Calculate Simpson's Diversity Index (1 - D)
  mutate(Diversity_Index = 1 - Simpsons_D)

pitfall_div <- merge(pitfall_div, overview_dataset, by = "Plot")

pitfall_div <- pitfall_div[pitfall_div$Diversity_Index != 0, ]

Model Evaluation:

Soil Carbon/Woodland Type Model - Non-Normal residuals

### Soil Carbon/Woodland Type Model - Normal residuals

###############################################################################################

# Soil Carbon/Woodland Type Model:

soil_carbon_model <- lme(fixed = Total_Carbon_Stock ~ Type * Layer_combined, 
                         random = ~ 1 | Site,  # Random intercepts and slopes for Type within Site
                         data = carbon_0_20, 
                         control = lmeControl(maxIter = 100, msMaxIter = 100, opt = "optim"))


# Soil Carbon Model R-squared:

soil_carbon_r_squared <- r.squaredGLMM(soil_carbon_model)

# Post-Hoc Results:

summary(soil_carbon_model)
## Linear mixed-effects model fit by REML
##   Data: carbon_0_20 
##        AIC      BIC    logLik
##   1278.591 1306.713 -629.2954
## 
## Random effects:
##  Formula: ~1 | Site
##         (Intercept) Residual
## StdDev:    14.93459 36.21726
## 
## Fixed effects:  Total_Carbon_Stock ~ Type * Layer_combined 
##                                                    Value Std.Error  DF
## (Intercept)                                     49.11022  12.54110 121
## TypeRegen                                      -13.15502  12.61891 121
## TypeUnforested (heath)                          -9.50266  13.54644 121
## TypeUnforested (remnant)                       -11.34849  12.80473 121
## Layer_combined0-20, O                            8.34979  13.09081 121
## TypeRegen:Layer_combined0-20, O                 11.39713  17.68497 121
## TypeUnforested (heath):Layer_combined0-20, O   -17.08628  19.10927 121
## TypeUnforested (remnant):Layer_combined0-20, O   3.81800  17.80741 121
##                                                  t-value p-value
## (Intercept)                                     3.915942  0.0001
## TypeRegen                                      -1.042485  0.2993
## TypeUnforested (heath)                         -0.701488  0.4843
## TypeUnforested (remnant)                       -0.886273  0.3772
## Layer_combined0-20, O                           0.637836  0.5248
## TypeRegen:Layer_combined0-20, O                 0.644453  0.5205
## TypeUnforested (heath):Layer_combined0-20, O   -0.894136  0.3730
## TypeUnforested (remnant):Layer_combined0-20, O  0.214405  0.8306
##  Correlation: 
##                                                (Intr) TypRgn TypU(h) TypU(r)
## TypeRegen                                      -0.520                       
## TypeUnforested (heath)                         -0.487  0.481                
## TypeUnforested (remnant)                       -0.511  0.507  0.473         
## Layer_combined0-20, O                          -0.507  0.499  0.469   0.489 
## TypeRegen:Layer_combined0-20, O                 0.372 -0.714 -0.344  -0.362 
## TypeUnforested (heath):Layer_combined0-20, O    0.346 -0.341 -0.710  -0.335 
## TypeUnforested (remnant):Layer_combined0-20, O  0.368 -0.365 -0.341  -0.719 
##                                                L_0-2O TR:L_O TypU(h):L_0-20,O
## TypeRegen                                                                    
## TypeUnforested (heath)                                                       
## TypeUnforested (remnant)                                                     
## Layer_combined0-20, O                                                        
## TypeRegen:Layer_combined0-20, O                -0.736                        
## TypeUnforested (heath):Layer_combined0-20, O   -0.683  0.503                 
## TypeUnforested (remnant):Layer_combined0-20, O -0.729  0.539  0.499          
## 
## Standardized Within-Group Residuals:
##        Min         Q1        Med         Q3        Max 
## -1.2472249 -0.6413349 -0.1636584  0.4138244  5.5070856 
## 
## Number of Observations: 131
## Number of Groups: 3
print(soil_carbon_r_squared)
##             R2m       R2c
## [1,] 0.05122386 0.1891091
###############################################################################################

# Soil Carbon/Woodland Type Model Residuals Normal Fit

residuals <- residuals(soil_carbon_model)

# Histogram to visualize residuals
hist(residuals, main = "Residuals Distribution", xlab = "Residuals", col = "lightblue", border = "black", breaks = 20)

# Q-Q plot to check for normality
qqnorm(residuals)
qqline(residuals)

# Shapiro-Wilk normality test
shapiro.test(residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals
## W = 0.80054, p-value = 4.619e-12
###############################################################################################

Seedling/Soil Carbon Model - Normal residuals

###############################################################################################


# Seedling/Soil Carbon Model:

upper10_seedling_model <- lme(fixed = sum_seedlings ~ Carbon_tonne_ha, 
                         random = ~ 1 | Site, 
                         data = regen_upper10)

upper10_est_seedling_model <- lme(fixed = sum_seedlings ~ Carbon_tonne_ha, 
                         random = ~ 1 | Site, 
                         data = regen_upper10_est)

# Seedling/Soil Carbon R-squared:

upper10_seedling_r_squared <- r.squaredGLMM(upper10_seedling_model)

upper10_est_seedling_r_squared <- r.squaredGLMM(upper10_est_seedling_model)


# Post-Hoc Results

summary(upper10_seedling_model)
## Linear mixed-effects model fit by REML
##   Data: regen_upper10 
##        AIC      BIC    logLik
##   209.4996 213.2773 -100.7498
## 
## Random effects:
##  Formula: ~1 | Site
##         (Intercept) Residual
## StdDev:    23.32199 32.94415
## 
## Fixed effects:  sum_seedlings ~ Carbon_tonne_ha 
##                    Value Std.Error DF   t-value p-value
## (Intercept)     90.73903 20.533222 17  4.419132  0.0004
## Carbon_tonne_ha -0.51144  0.420827 17 -1.215323  0.2409
##  Correlation: 
##                 (Intr)
## Carbon_tonne_ha -0.669
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -1.99297649 -0.52684978  0.04632305  0.63364651  1.38446672 
## 
## Number of Observations: 21
## Number of Groups: 3
print(upper10_seedling_r_squared)
##            R2m      R2c
## [1,] 0.0578357 0.372375
summary(upper10_est_seedling_model)
## Linear mixed-effects model fit by REML
##   Data: regen_upper10_est 
##        AIC     BIC    logLik
##   172.6892 176.022 -82.34459
## 
## Random effects:
##  Formula: ~1 | Site
##         (Intercept) Residual
## StdDev:    26.43445 18.74382
## 
## Fixed effects:  sum_seedlings ~ Carbon_tonne_ha 
##                   Value Std.Error DF    t-value p-value
## (Intercept)     41.4879 17.738754 15  2.3388286  0.0336
## Carbon_tonne_ha -0.0868  0.246273 15 -0.3524518  0.7294
##  Correlation: 
##                 (Intr)
## Carbon_tonne_ha -0.447
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -2.69250319 -0.43396292 -0.03079282  0.60978544  1.48404159 
## 
## Number of Observations: 19
## Number of Groups: 3
print(upper10_est_seedling_r_squared)
##              R2m      R2c
## [1,] 0.002990246 0.666435
###############################################################################################

# Seedling/Soil Carbon Model Residuals Normal Fit

residuals_upper10_seedlings <- residuals(upper10_seedling_model)

residuals_upper10_est_seedlings <- residuals(upper10_est_seedling_model)


# Histogram to visualize residuals
hist(residuals_upper10_seedlings, main = "Residuals Distribution", xlab = "Residuals", col = "lightblue", border = "black")

# Histogram to visualize residuals
hist(residuals_upper10_est_seedlings, main = "Residuals Distribution", xlab = "Residuals", col = "lightblue", border = "black")

# Q-Q plot to check for normality
qqnorm(residuals_upper10_seedlings)
qqline(residuals_upper10_seedlings)

# Shapiro-Wilk normality test
shapiro.test(residuals_upper10_seedlings)
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals_upper10_seedlings
## W = 0.96306, p-value = 0.5797
# Q-Q plot to check for normality
qqnorm(residuals_upper10_est_seedlings)
qqline(residuals_upper10_est_seedlings)

# Shapiro-Wilk normality test
shapiro.test(residuals_upper10_est_seedlings)
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals_upper10_est_seedlings
## W = 0.91636, p-value = 0.09694
###############################################################################################

Plant Diversity Model - Non-Normal residuals

###############################################################################################

# Plant Diversity Model

diversity_model <- lme(fixed = Diversity_Index ~ Type, 
                         random = ~ 1 | Site, 
                         data = plant_div,
                       control = lmeControl(maxIter = 100, msMaxIter = 100, opt = "optim"))


# Plant Diversity Model R-squared:

plant_diversity_r_squared <- r.squaredGLMM(diversity_model)


summary(diversity_model)
## Linear mixed-effects model fit by REML
##   Data: plant_div 
##         AIC       BIC   logLik
##   -79.72206 -65.89767 45.86103
## 
## Random effects:
##  Formula: ~1 | Site
##         (Intercept)  Residual
## StdDev: 0.003486986 0.1201534
## 
## Fixed effects:  Diversity_Index ~ Type 
##                               Value  Std.Error DF   t-value p-value
## (Intercept)               0.6861663 0.02629680 72 26.093143  0.0000
## TypeRegen                -0.0356075 0.03708015 72 -0.960284  0.3401
## TypeUnforested (heath)    0.0490610 0.04061927 72  1.207826  0.2311
## TypeUnforested (remnant)  0.0139425 0.03708015 72  0.376009  0.7080
##  Correlation: 
##                          (Intr) TypRgn TypU(h)
## TypeRegen                -0.705               
## TypeUnforested (heath)   -0.644  0.456        
## TypeUnforested (remnant) -0.705  0.500  0.456 
## 
## Standardized Within-Group Residuals:
##        Min         Q1        Med         Q3        Max 
## -3.8956063 -0.5227355  0.1062503  0.6021995  1.6875788 
## 
## Number of Observations: 78
## Number of Groups: 3
print(plant_diversity_r_squared)
##             R2m        R2c
## [1,] 0.05584031 0.05663483
###############################################################################################

# Plant Diversity Model Residuals Normal Fit

residuals_plant_diversity <- residuals(diversity_model)

# Histogram to visualize residuals
hist(residuals_plant_diversity, main = "Residuals Distribution", xlab = "Residuals", col = "lightblue", border = "black")

# Q-Q plot to check for normality
qqnorm(residuals_plant_diversity)
qqline(residuals_plant_diversity)

# Shapiro-Wilk normality test
shapiro.test(residuals_plant_diversity)
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals_plant_diversity
## W = 0.94101, p-value = 0.001287
###############################################################################################

Plant Richness Model - Not normal residuals

# Mixed model (if it fits)
richness_model_lme <- glmer(Richness ~ Type + (1 | Site), data = plant_rich, family = poisson)
## boundary (singular) fit: see help('isSingular')
# Fixed model
richness_model_lm <- lm(Richness ~ Type, data = plant_rich)

# Compare (will likely show negligible gain in model fit)
anova(richness_model_lme, richness_model_lm)
## Data: plant_rich
## Models:
## richness_model_lme: Richness ~ Type + (1 | Site)
## richness_model_lm: Richness ~ Type
##                    npar    AIC    BIC  logLik -2*log(L) Chisq Df Pr(>Chisq)
## richness_model_lme    5 383.76 395.54 -186.88    373.76                    
## richness_model_lm     5 378.76 390.54 -184.38    368.76 5.002  0
richness_model_lm <- glm(Richness ~ Type, data = plant_rich, family = poisson())

deviance <- deviance(richness_model_lm)
df_resid <- df.residual(richness_model_lm)
dispersion_ratio <- deviance / df_resid
dispersion_ratio
## [1] 0.6040578
###############################################################################################

# Plant Richness Model

richness_model <- glm(Richness ~ Type, data = plant_rich, family = poisson())

# Plant Diversity Model R-squared:

plant_richness_r_squared <- r.squaredGLMM(richness_model)
## Warning: the null model is only correct if all the variables it uses are identical 
## to those used in fitting the original model.
summary(richness_model)
## 
## Call:
## glm(formula = Richness ~ Type, family = poisson(), data = plant_rich)
## 
## Coefficients:
##                          Estimate Std. Error z value Pr(>|z|)    
## (Intercept)              2.339973   0.067729  34.549   <2e-16 ***
## TypeRegen                0.168465   0.091995   1.831   0.0671 .  
## TypeUnforested (heath)   0.014572   0.104481   0.139   0.8891    
## TypeUnforested (remnant) 0.009132   0.095565   0.096   0.9239    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 49.301  on 77  degrees of freedom
## Residual deviance: 44.700  on 74  degrees of freedom
## AIC: 381.76
## 
## Number of Fisher Scoring iterations: 4
print(plant_richness_r_squared)
##                  R2m        R2c
## delta     0.05395323 0.05395323
## lognormal 0.05624468 0.05624468
## trigamma  0.05165377 0.05165377
###############################################################################################

# Plant Diversity Model Residuals Normal Fit

residuals_plant_richness <- residuals(richness_model)

# Histogram to visualize residuals
hist(residuals_plant_richness, main = "Residuals Distribution", xlab = "Residuals", col = "lightblue", border = "black")

# Q-Q plot to check for normality
qqnorm(residuals_plant_richness)
qqline(residuals_plant_richness)

# Shapiro-Wilk normality test
shapiro.test(residuals_plant_richness)
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals_plant_richness
## W = 0.93667, p-value = 0.0007625
###############################################################################################

Soil Invertebrate Diversity Model - Normal Residuals

###############################################################################################

# Soil Invertebrate Diversity Model

soil_invert_diversity_model <- lme(fixed = Diversity_Index ~ Type, 
                         random = ~ 1 | Site, 
                         data = soil_invert_div,
                         , 
                         control = lmeControl(maxIter = 100, msMaxIter = 100, opt = "optim"))

# Soil Invertebrate Diversity Model R-squared:

soil_invert_diversity_r_squared <- r.squaredGLMM(soil_invert_diversity_model)


summary(soil_invert_diversity_model)
## Linear mixed-effects model fit by REML
##   Data: soil_invert_div 
##         AIC      BIC   logLik
##   -3.784492 7.187356 7.892246
## 
## Random effects:
##  Formula: ~1 | Site
##         (Intercept)  Residual
## StdDev:  0.00152749 0.1826833
## 
## Fixed effects:  Diversity_Index ~ Type 
##                               Value  Std.Error DF   t-value p-value
## (Intercept)               0.3659797 0.04883611 45  7.494038  0.0000
## TypeRegen                 0.0949750 0.07036316 45  1.349784  0.1838
## TypeUnforested (heath)   -0.0654500 0.07563807 45 -0.865305  0.3915
## TypeUnforested (remnant) -0.0519503 0.07036316 45 -0.738317  0.4642
##  Correlation: 
##                          (Intr) TypRgn TypU(h)
## TypeRegen                -0.694               
## TypeUnforested (heath)   -0.645  0.448        
## TypeUnforested (remnant) -0.694  0.481  0.448 
## 
## Standardized Within-Group Residuals:
##        Min         Q1        Med         Q3        Max 
## -2.5234455 -0.6473544  0.1723194  0.5872532  1.9679246 
## 
## Number of Observations: 50
## Number of Groups: 2
print(soil_invert_diversity_r_squared)
##           R2m       R2c
## [1,] 0.106536 0.1065985
###############################################################################################

# Soil Invertebrate Diversity Model Residuals Normal Fit

residuals_soil_invert_diversity <- residuals(soil_invert_diversity_model)

# Histogram to visualize residuals
hist(residuals_soil_invert_diversity, main = "Residuals Distribution", xlab = "Residuals", col = "lightblue", border = "black")

# Q-Q plot to check for normality
qqnorm(residuals_soil_invert_diversity)
qqline(residuals_soil_invert_diversity)

# Shapiro-Wilk normality test
shapiro.test(residuals_soil_invert_diversity)
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals_soil_invert_diversity
## W = 0.98845, p-value = 0.9036
###############################################################################################

Soil Invertebrate Richness Model - Not normal residuals

# Mixed model (if it fits)
soil_richness_invert_model_lme <- glmer(Richness ~ Type + (1 | Site), data = soil_invert_rich, family = poisson)
## boundary (singular) fit: see help('isSingular')
soil_invert_richness_model_lm <- glm(Richness ~ Type, 
                                      data = soil_invert_rich, 
                                      family = poisson)

# Compare (will likely show negligible gain in model fit)
anova(soil_richness_invert_model_lme, soil_invert_richness_model_lm)
## Data: soil_invert_rich
## Models:
## soil_invert_richness_model_lm: Richness ~ Type
## soil_richness_invert_model_lme: Richness ~ Type + (1 | Site)
##                                npar    AIC    BIC logLik -2*log(L) Chisq Df
## soil_invert_richness_model_lm     4 168.72 176.37 -80.36    160.72         
## soil_richness_invert_model_lme    5 170.72 180.28 -80.36    160.72     0  1
##                                Pr(>Chisq)
## soil_invert_richness_model_lm            
## soil_richness_invert_model_lme          1
deviance <- deviance(soil_invert_richness_model_lm)
df_resid <- df.residual(soil_invert_richness_model_lm)
dispersion_ratio <- deviance / df_resid
dispersion_ratio
## [1] 0.3068677
###############################################################################################

# Soil Invertebrate Richness Model

soil_invert_richness_model <- glm(Richness ~ Type, 
                                      data = soil_invert_rich, 
                                      family = poisson)

# Soil Invertebrate Richness Model R-squared:

soil_invert_richness_r_squared <- r.squaredGLMM(soil_invert_richness_model)
## Warning: the null model is only correct if all the variables it uses are identical 
## to those used in fitting the original model.
summary(soil_invert_richness_model)
## 
## Call:
## glm(formula = Richness ~ Type, family = poisson, data = soil_invert_rich)
## 
## Coefficients:
##                          Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                1.3312     0.1374   9.692   <2e-16 ***
## TypeRegen                 -0.2853     0.2142  -1.332   0.1830    
## TypeUnforested (heath)    -0.3380     0.2364  -1.429   0.1529    
## TypeUnforested (remnant)  -0.3997     0.2217  -1.802   0.0715 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 18.163  on 49  degrees of freedom
## Residual deviance: 14.116  on 46  degrees of freedom
## AIC: 168.72
## 
## Number of Fisher Scoring iterations: 4
print(soil_invert_richness_r_squared)
##                  R2m        R2c
## delta     0.07155917 0.07155917
## lognormal 0.08198368 0.08198368
## trigamma  0.06107931 0.06107931
###############################################################################################

# Soil Invertebrate Richness Model Residuals Normal Fit

residuals_soil_invert_richness <- residuals(soil_invert_richness_model)

# Histogram to visualize residuals
hist(residuals_soil_invert_richness, main = "Residuals Distribution", xlab = "Residuals", col = "lightblue", border = "black")

# Q-Q plot to check for normality
qqnorm(residuals_soil_invert_richness)
qqline(residuals_soil_invert_richness)

# Shapiro-Wilk normality test
shapiro.test(residuals_soil_invert_richness)
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals_soil_invert_richness
## W = 0.94179, p-value = 0.01587
###############################################################################################

Pitfall Invertebrate Diversity Model - Not Normal residuals

###############################################################################################

# Pitfall Invertebrate Diversity Model

pitfall_diversity_model <- lme(fixed = Diversity_Index ~ Type, 
                         random = ~ 1 | Site, 
                         data = pitfall_div,
                         , 
                         control = lmeControl(maxIter = 100, msMaxIter = 100, opt = "optim"))



# Soil Invertebrate Diversity Model R-squared:

pitfall_invert_diversity_r_squared <- r.squaredGLMM(pitfall_diversity_model)


summary(pitfall_diversity_model)
## Linear mixed-effects model fit by REML
##   Data: pitfall_div 
##        AIC      BIC   logLik
##   -6.73401 6.756961 9.367005
## 
## Random effects:
##  Formula: ~1 | Site
##         (Intercept)  Residual
## StdDev:  0.05952625 0.1914478
## 
## Fixed effects:  Diversity_Index ~ Type 
##                               Value  Std.Error DF   t-value p-value
## (Intercept)               0.6610537 0.05491415 68 12.037949  0.0000
## TypeRegen                -0.0025667 0.06135071 68 -0.041836  0.9668
## TypeUnforested (heath)    0.0776850 0.06671556 68  1.164421  0.2483
## TypeUnforested (remnant) -0.0045245 0.05983130 68 -0.075620  0.9399
##  Correlation: 
##                          (Intr) TypRgn TypU(h)
## TypeRegen                -0.543               
## TypeUnforested (heath)   -0.500  0.448        
## TypeUnforested (remnant) -0.558  0.499  0.459 
## 
## Standardized Within-Group Residuals:
##        Min         Q1        Med         Q3        Max 
## -2.8213320 -0.4327584  0.2168915  0.6971038  1.6957427 
## 
## Number of Observations: 74
## Number of Groups: 3
print(pitfall_invert_diversity_r_squared)
##             R2m       R2c
## [1,] 0.02427668 0.1102898
###############################################################################################

# Soil Invertebrate Diversity Model Residuals Normal Fit

residuals_pitfall_invert_diversity <- residuals(pitfall_diversity_model)

# Histogram to visualize residuals
hist(residuals_pitfall_invert_diversity, main = "Residuals Distribution", xlab = "Residuals", col = "lightblue", border = "black")

# Q-Q plot to check for normality
qqnorm(residuals_pitfall_invert_diversity)
qqline(residuals_pitfall_invert_diversity)

# Shapiro-Wilk normality test
shapiro.test(residuals_pitfall_invert_diversity)
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals_pitfall_invert_diversity
## W = 0.93568, p-value = 0.0009668
###############################################################################################

Pitfall Invertebrate Richness Model - Normal Residuals

###############################################################################################

# Pitfall Invertebrate Richness Model

pitfall_richness_model <- glmer(Richness ~ Type + 
                         (1 | Site),
                         family = poisson(),
                         data = pitfall_rich)

# Soil Invertebrate Richness Model R-squared:

pitfall_invert_richness_r_squared <- r.squaredGLMM(pitfall_richness_model)
## Warning: the null model is only correct if all the variables it uses are identical 
## to those used in fitting the original model.
summary(pitfall_richness_model)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: Richness ~ Type + (1 | Site)
##    Data: pitfall_rich
## 
##       AIC       BIC    logLik -2*log(L)  df.resid 
##     340.0     351.5    -165.0     330.0        69 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.78501 -0.63764 -0.02627  0.67942  2.50383 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Site   (Intercept) 0.02481  0.1575  
## Number of obs: 74, groups:  Site, 3
## 
## Fixed effects:
##                          Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                1.7718     0.1278  13.861   <2e-16 ***
## TypeRegen                 -0.2433     0.1348  -1.806    0.071 .  
## TypeUnforested (heath)    -0.1389     0.1540  -0.902    0.367    
## TypeUnforested (remnant)  -0.1735     0.1337  -1.298    0.194    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) TypRgn TypU(h)
## TypeRegen   -0.463               
## TypUnf(hth) -0.397  0.385        
## TypU(rmnnt) -0.465  0.443  0.391
print(pitfall_invert_richness_r_squared)
##                  R2m       R2c
## delta     0.03992091 0.1497831
## lognormal 0.04304071 0.1614886
## trigamma  0.03673198 0.1378183
###############################################################################################

# Soil Invertebrate Richness Model Residuals Normal Fit

residuals_pitfall_invert_richness <- residuals(pitfall_richness_model)

# Histogram to visualize residuals
hist(residuals_pitfall_invert_richness, main = "Residuals Distribution", xlab = "Residuals", col = "lightblue", border = "black")

# Q-Q plot to check for normality
qqnorm(residuals_pitfall_invert_richness)
qqline(residuals_pitfall_invert_richness)

# Shapiro-Wilk normality test
shapiro.test(residuals_pitfall_invert_richness)
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals_pitfall_invert_richness
## W = 0.96675, p-value = 0.04836
###############################################################################################

Analysis:

Soil Carbon:

Extracting model predicted estimates:

library(nlme)

# Subset the data for Organic Soil (Layer_combined == '0-20, O')
organic_data <- carbon_0_20[carbon_0_20$Layer_combined == "0-20, O", ]

# Subset the data for Mineral Soil (Layer_combined == '0-10, M')
mineral_data <- carbon_0_20[carbon_0_20$Layer_combined == "0-10, M", ]

library(nlme)

organic_model <- lme(
  Total_Carbon_Stock ~ Type,
  random = ~ 1 | Site,
  data = organic_data
)

mineral_model <- lme(
  Total_Carbon_Stock ~ Type,
  random = ~ 1 | Site,
  data = mineral_data
)


organic_confints <- intervals(organic_model, level = 0.95)$fixed
mineral_confints <- intervals(mineral_model, which = "fixed", level = 0.95)

mineral_fixed <- mineral_confints$fixed

# Step 1: Get the intercept from the organic and mineral models
organic_intercept <- organic_confints["(Intercept)", "est."]
mineral_intercept <- mineral_fixed["(Intercept)", "est."]  # or "Estimate" depending on version

# Step 2: Add the intercept to other treatment-level estimates
organic_adjusted <- organic_confints
organic_adjusted[-1, ] <- organic_adjusted[-1, ] + organic_intercept  # Skip the intercept row

# Adjust the rest
mineral_adjusted <- mineral_fixed
mineral_adjusted[rownames(mineral_adjusted) != "(Intercept)", ] <-
  mineral_adjusted[rownames(mineral_adjusted) != "(Intercept)", ] + mineral_intercept

mineral_adjusted
##                             lower     est.    upper
## (Intercept)              36.97371 45.69073 54.40775
## TypeRegen                21.53700 33.68209 45.82717
## TypeUnforested (heath)   26.08765 39.10718 52.12671
## TypeUnforested (remnant) 22.01451 34.34224 46.66997
## attr(,"label")
## [1] "Fixed effects:"
# Step 3: Rebuild the cleaned-up dataframes

organic_df <- data.frame(
  Type = c("Mature", "Regen", "Unforested (heath)", "Unforested (remnant)"),
  Soil_Type = "Organic",
  Lower = organic_adjusted[, "lower"],
  Estimate = organic_adjusted[, "est."],
  Upper = organic_adjusted[, "upper"]
)

mineral_df <- data.frame(
  Type = c("Mature", "Regen", "Unforested (heath)", "Unforested (remnant)"),
  Soil_Type = "Mineral",
  Lower = mineral_adjusted[, "lower"],
  Estimate = mineral_adjusted[, "est."],
  Upper = mineral_adjusted[, "upper"]
)

# Combine
combined_df <- rbind(organic_df, mineral_df)
library(dplyr)

# Sort just to be safe
stacked_df <- combined_df %>%
  arrange(Type, Soil_Type) %>%
  mutate(
    Soil_Type = factor(Soil_Type, levels = c("Organic", "Mineral"))
  ) %>%
  group_by(Type) %>%
  arrange(Soil_Type, .by_group = TRUE) %>%
  mutate(
    ymin = ifelse(Soil_Type == "Organic", 0, lag(Estimate)),
    ymax = ymin + Estimate,
    
    # For error bars
    error_lower = ifelse(Soil_Type == "Organic", Lower, lag(Estimate) + Lower),
    error_upper = ifelse(Soil_Type == "Organic", Upper, lag(Estimate) + Upper)
  )


# Reorder the Soil_Type factor levels
stacked_df$Soil_Type <- factor(stacked_df$Soil_Type, levels = c("Mineral", "Organic"))  # Change order here

stacked_df
## # A tibble: 8 × 9
## # Groups:   Type [4]
##   Type        Soil_Type Lower Estimate Upper  ymin  ymax error_lower error_upper
##   <chr>       <fct>     <dbl>    <dbl> <dbl> <dbl> <dbl>       <dbl>       <dbl>
## 1 Mature      Organic   21.7      56.7  91.7   0    56.7       21.7         91.7
## 2 Mature      Mineral   37.0      45.7  54.4  56.7 102.        93.7        111. 
## 3 Regen       Organic   24.0      55.6  87.1   0    55.6       24.0         87.1
## 4 Regen       Mineral   21.5      33.7  45.8  55.6  89.2       77.1        101. 
## 5 Unforested… Organic   -3.72     30.6  64.9   0    30.6       -3.72        64.9
## 6 Unforested… Mineral   26.1      39.1  52.1  30.6  69.7       56.7         82.7
## 7 Unforested… Organic   18.2      49.7  81.2   0    49.7       18.2         81.2
## 8 Unforested… Mineral   22.0      34.3  46.7  49.7  84.1       71.7         96.4
organic_estimates <- stacked_df %>%
  dplyr::filter(Soil_Type == "Organic") %>%  # Filter only organic estimates
  dplyr::select(Type, Estimate)  # Select Type and Estimate columns

# View the organic estimates
print(organic_estimates)
## # A tibble: 4 × 2
## # Groups:   Type [4]
##   Type                 Estimate
##   <chr>                   <dbl>
## 1 Mature                   56.7
## 2 Regen                    55.6
## 3 Unforested (heath)       30.6
## 4 Unforested (remnant)     49.7
# Join the organic estimates with the raw carbon data
carbon_0_20_adjusted <- carbon_0_20 %>%
  left_join(organic_estimates, by = "Type")  # Join by 'Type' (treatment)

# View the adjusted data to confirm the join
head(carbon_0_20_adjusted)
## # A tibble: 6 × 6
## # Groups:   Plot, Type, Layer_combined [6]
##    Plot Type                 Layer_combined Site     Total_Carbon_Stock Estimate
##   <dbl> <chr>                <chr>          <chr>                 <dbl>    <dbl>
## 1     1 Unforested (remnant) 0-20, O        Glen Ma…               24.5     49.7
## 2     2 Unforested (remnant) 0-20, O        Glen Ma…               31.2     49.7
## 3     3 Mature               0-20, O        Glen Ma…              193.      56.7
## 4     4 Mature               0-10, M        Glen Ma…               85.0     56.7
## 5     5 Unforested (remnant) 0-10, M        Glen Ma…               52.3     49.7
## 6     5 Unforested (remnant) 0-20, O        Glen Ma…               72.3     49.7
# Adjust the mineral carbon stock values by adding the organic estimate
carbon_0_20_adjusted <- carbon_0_20_adjusted %>%
  mutate(
    Total_Carbon_Stock_adjusted = ifelse(
      Layer_combined == 'Mineral',  # If it's mineral soil
      Total_Carbon_Stock + Estimate,  # Add the organic estimate
      Total_Carbon_Stock  # Keep the original carbon stock for organic soil
    )
  )

# View the adjusted data
head(carbon_0_20_adjusted)
## # A tibble: 6 × 7
## # Groups:   Plot, Type, Layer_combined [6]
##    Plot Type                 Layer_combined Site     Total_Carbon_Stock Estimate
##   <dbl> <chr>                <chr>          <chr>                 <dbl>    <dbl>
## 1     1 Unforested (remnant) 0-20, O        Glen Ma…               24.5     49.7
## 2     2 Unforested (remnant) 0-20, O        Glen Ma…               31.2     49.7
## 3     3 Mature               0-20, O        Glen Ma…              193.      56.7
## 4     4 Mature               0-10, M        Glen Ma…               85.0     56.7
## 5     5 Unforested (remnant) 0-10, M        Glen Ma…               52.3     49.7
## 6     5 Unforested (remnant) 0-20, O        Glen Ma…               72.3     49.7
## # ℹ 1 more variable: Total_Carbon_Stock_adjusted <dbl>
# Recode the 'Layer_combined' column to use 'Mineral' and 'Organic'
carbon_0_20_adjusted <- carbon_0_20_adjusted %>%
  mutate(
    Layer_combined = recode(Layer_combined, 
                            '0-10, M' = 'Mineral',  # Change '0-10, M' to 'Mineral'
                            '0-20, O' = 'Organic')   # Optionally, recode 'Organic' if needed
  )

Plotting soil carbon stock

stacked_df$Treatment_Layer <- paste(stacked_df$Type, stacked_df$Soil_Type, sep = " / ")
carbon_0_20_adjusted$Treatment_Layer <- paste(carbon_0_20_adjusted$Type, carbon_0_20_adjusted$Layer_combined, sep = " / ")


library(colorspace)
palette.colors()
## [1] "#000000" "#E69F00" "#56B4E9" "#009E73" "#F0E442" "#0072B2" "#D55E00"
## [8] "#CC79A7" "#999999"
# Base woodland treatment colors
woodland_base <- c(
  "Mature" = "#D55E00", 
  "Regen" = "#009E73", 
  "Unforested (heath)" = "#56B4E9", 
  "Unforested (remnant)" = "#CC79A7"
)

# Build combined color palette
layer_colors <- c()
for (type in names(woodland_base)) {
  layer_colors[paste0(type, " / Organic")] <- woodland_base[type]
  layer_colors[paste0(type, " / Mineral")] <- darken(woodland_base[type], amount = 0.2)
}
library(ggplot2)

soil_woodland_graph <- ggplot(stacked_df, aes(x = Type, y = Estimate, fill = Treatment_Layer)) +
  geom_bar(stat = "identity", position = "stack", color = "black", size = 0.5, width = 0.6, alpha = 0.7) +
  geom_errorbar(aes(ymin = error_lower, ymax = error_upper),
                position = "identity", width = 0.3, color = "black", size = 0.5) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "black") +
  scale_y_reverse() +
  scale_shape_manual(values = c("Organic" = 21, "Mineral" = 24)) +
  scale_fill_manual(values = layer_colors) +
  guides(
    fill = guide_legend(override.aes = list(shape = c(21, 24)), reverse = TRUE),
    shape = "none"  # Hide the separate shape legend
  ) +
  theme_bw() +
  theme(
  axis.text.x = element_text(angle = 45, hjust = 1),  # <-- added line
  strip.background = element_rect(fill = "grey", color = "black", size = 1),
  strip.text = element_text(color = "black", size = 12),
  panel.spacing = unit(1, "lines"),
  panel.background = element_rect(fill = "white", color = "black"),
  plot.margin = margin(10, 10, 10, 10)
) +
  labs(
    x = "Treatment Type",
    y = "Carbon Stock (tonne/ha)",
    fill = "Forest Treatment / Soil Layer",
    title = "Stacked Organic & Mineral Soil Carbon by Treatment"
  ) +
  geom_jitter(data = carbon_0_20_adjusted, 
              aes(x = Type,
                  y = ifelse(Layer_combined == 'Mineral', 
                             Total_Carbon_Stock_adjusted, 
                             Total_Carbon_Stock),
                  fill = Treatment_Layer,
                  shape = Layer_combined),
              width = 0.15, 
              height = 0, 
              size = 1.5, 
              color = "black", 
              stroke = 0.8, 
              alpha = 0.5)
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## Warning: The `size` argument of `element_rect()` is deprecated as of ggplot2 3.4.0.
## ℹ Please use the `linewidth` argument instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
soil_woodland_graph

save_plot(soil_woodland_graph, "soil_carbon_graph", plot_title = "Carbon stock by Forest Treatment")
## png 
##   2

Plotting Organic Depth Horizon

# Modify the O_layer_depth column
depth_dataset <- depth_dataset %>%
  mutate(
    O_layer_depth = case_when(
      O_layer_depth == ">20" ~ 25,        # Replace ">20" with 25
      O_layer_depth == "(=20)" ~ 20,      # Replace "(=20)" with 20
      TRUE ~ as.numeric(O_layer_depth)   # Convert remaining values to numeric
    )
  ) %>%
  filter(!is.na(O_layer_depth))  # Remove rows where O_layer_depth is NA
## Warning: There was 1 warning in `mutate()`.
## ℹ In argument: `O_layer_depth = case_when(...)`.
## Caused by warning:
## ! NAs introduced by coercion
depth_dataset <- depth_dataset[, c("Plot", "O_layer_depth", "Site")]

depth_dataset <- merge(depth_dataset, overview_dataset, by = c("Plot", "Site"))

depth_dataset <- depth_dataset[, c("Plot", "O_layer_depth", "Type", "Site")]

average_depth_results <- depth_dataset %>%
  group_by(Type, Site) %>%
  summarise(
    Average_O_layer_depth = mean(O_layer_depth, na.rm = TRUE),  # Mean depth
    SE_O_layer_depth = sd(O_layer_depth, na.rm = TRUE) / sqrt(n()),  # Standard Error
    Lower_CI_depth = Average_O_layer_depth - 1.96 * SE_O_layer_depth,  # Lower 95% CI
    Upper_CI_depth = Average_O_layer_depth + 1.96 * SE_O_layer_depth   # Upper 95% CI
  )
## `summarise()` has grouped output by 'Type'. You can override using the
## `.groups` argument.
average_plot_depth <- depth_dataset %>%
  group_by(Plot, Type, Site) %>%
  summarise(
    Plot_O_layer_depth = mean(O_layer_depth, na.rm = TRUE))
## `summarise()` has grouped output by 'Plot', 'Type'. You can override using the
## `.groups` argument.
woodland_base <- c(
  "Mature" = "#D55E00", 
  "Regen" = "#009E73", 
  "Unforested (heath)" = "#56B4E9", 
  "Unforested (remnant)" = "#CC79A7"
)

depth_graph <- ggplot(average_depth_results, aes(x = Type, y = Average_O_layer_depth, fill = Type)) + 
  geom_bar(stat = "identity", position = "stack", color = "black", size = 0.5, width = 0.6, alpha = 0.7) +
  scale_fill_manual(name = "Forest Treatment",
                    values = c(
  "Mature" = "#D55E00", 
  "Regen" = "#009E73", 
  "Unforested (heath)" = "#56B4E9", 
  "Unforested (remnant)" = "#CC79A7"
  )) +
  geom_errorbar(data = average_depth_results, 
                aes(ymin = Lower_CI_depth, ymax = Upper_CI_depth, x = Type), 
                width = 0.2, color = "black", size = 0.5) + 
  geom_hline(yintercept = 0, linetype = "dashed", color = "black") +
  scale_y_reverse() +
  theme_minimal() +  # Use a minimal theme with a white background
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1),  # Rotate x-axis labels
    strip.background = element_rect(fill = "grey", color = "black", size = 1),  # Grey background for facet labels
    strip.text = element_text(color = "black", size = 12),  # Text color inside the grey box
    panel.spacing = unit(1, "lines"),  # Add spacing between panels if needed
    panel.background = element_rect(fill = "white", color = "black"),  # Apply white background with black border to panels
    plot.margin = margin(10, 10, 10, 10) # Add margins around the plot if necessary
  ) +# Reverse y-axis for better visibility
  labs(
    x = "Forest Treatment",  # Custom x-axis label
    y = "Organic Layer Depth (cm)",  # Custom y-axis label
    fill = "Soil Layer",  # Custom legend title
    title = "Organic Horizon Depth across Sites"  # Overall title
  ) +  # Create facets based on 'Site'
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1)) +
  facet_wrap(~Site)

depth_graph

save_plot(depth_graph, "depth_graph")
## png 
##   2
# Perform one-way ANOVA for each site
anova_results <- depth_dataset %>%
  do({
    model <- aov(O_layer_depth ~ Type, data = .)  # Perform ANOVA within each site
    tidy(model)  # Use the tidy() function from broom to get results in a nice format
  })

# View the results of the ANOVA for each site
anova_results
## # A tibble: 2 × 6
##   term         df  sumsq meansq statistic p.value
##   <chr>     <dbl>  <dbl>  <dbl>     <dbl>   <dbl>
## 1 Type          3   429.  143.       1.58   0.195
## 2 Residuals   351 31836.   90.7     NA     NA
# Perform Tukey's HSD test for each site where ANOVA was significant
tukey_results <- depth_dataset %>%
  group_by(Site) %>%
  do({
    model <- aov(O_layer_depth ~ Type, data = .)  # Fit the ANOVA model
    tukey <- TukeyHSD(model)  # Perform Tukey's HSD post-hoc test
    tidy(tukey)  # Get results in a tidy format
  })

# View Tukey HSD results for each site
tukey_results
## # A tibble: 18 × 8
## # Groups:   Site [3]
##    Site        term  contrast null.value estimate conf.low conf.high adj.p.value
##    <chr>       <chr> <chr>         <dbl>    <dbl>    <dbl>     <dbl>       <dbl>
##  1 Glen Cluan… Type  Regen-M…          0  2.26      -3.01     7.53        0.680 
##  2 Glen Cluan… Type  Unfores…          0  1.22      -4.43     6.87        0.943 
##  3 Glen Cluan… Type  Unfores…          0  3.23      -2.00     8.46        0.377 
##  4 Glen Cluan… Type  Unfores…          0 -1.04      -6.53     4.45        0.960 
##  5 Glen Cluan… Type  Unfores…          0  0.974     -4.08     6.03        0.958 
##  6 Glen Cluan… Type  Unfores…          0  2.01      -3.44     7.47        0.771 
##  7 Glen Loyne  Type  Regen-M…          0  1.26      -3.33     5.86        0.891 
##  8 Glen Loyne  Type  Unfores…          0  1.93      -3.05     6.91        0.745 
##  9 Glen Loyne  Type  Unfores…          0  1.25      -3.31     5.82        0.891 
## 10 Glen Loyne  Type  Unfores…          0  0.666     -4.25     5.58        0.985 
## 11 Glen Loyne  Type  Unfores…          0 -0.00840   -4.50     4.48        1.00  
## 12 Glen Loyne  Type  Unfores…          0 -0.674     -5.56     4.21        0.984 
## 13 Glen Mallie Type  Regen-M…          0  0.866     -5.13     6.86        0.982 
## 14 Glen Mallie Type  Unfores…          0 -6.5      -13.7      0.657       0.0891
## 15 Glen Mallie Type  Unfores…          0  0.883     -5.01     6.77        0.980 
## 16 Glen Mallie Type  Unfores…          0 -7.37     -14.6     -0.0871      0.0462
## 17 Glen Mallie Type  Unfores…          0  0.0167    -6.02     6.06        1.00  
## 18 Glen Mallie Type  Unfores…          0  7.38       0.188   14.6         0.0420

Seedling Establishmet on soil carbon:

Extracting and plotting model predicted data:

regen_upper10_graph    <- make_seedling_plot(regen_upper10, "Regenerating", color_line = "#009E73", color_point = "black", custom_title = "Total Seedlings vs Topsoil SOC")



regen_upper10_est_graph    <- make_seedling_plot(regen_upper10_est, "Regenerating", color_line = "#009E73", color_point = "black", custom_title = "Established Seedlings vs Topsoil SOC")



combined_regen_graph <- regen_upper10_graph / regen_upper10_est_graph

combined_regen_graph

save_plot(combined_regen_graph, "seedling_regen_graph")
## png 
##   2

Plant Diversity:

Extracting model predicted data

predicted_plant_div_with_ci <- intervals(diversity_model, level = 0.95)

# Extract just the fixed effects into a data frame
plant_div_ci_df <- as.data.frame(predicted_plant_div_with_ci$fixed)

# Optionally, name the columns
colnames(plant_div_ci_df) <- c("CI_lower", "Estimate", "CI_upper")

# Add a column for the effect names (e.g., Intercept, Carbon, etc.)
plant_div_ci_df$Effect <- rownames(predicted_plant_div_with_ci$fixed)

# Reorder for readability
plant_div_ci_df <- plant_div_ci_df[, c("Effect", "Estimate", "CI_lower", "CI_upper")]

# View it
print(plant_div_ci_df)
##                                            Effect    Estimate    CI_lower
## (Intercept)                           (Intercept)  0.68616626  0.63374454
## TypeRegen                               TypeRegen -0.03560748 -0.10952542
## TypeUnforested (heath)     TypeUnforested (heath)  0.04906101 -0.03191202
## TypeUnforested (remnant) TypeUnforested (remnant)  0.01394246 -0.05997548
##                            CI_upper
## (Intercept)              0.73858798
## TypeRegen                0.03831045
## TypeUnforested (heath)   0.13003405
## TypeUnforested (remnant) 0.08786039
# Extract the intercept value
intercept_estimate <- plant_div_ci_df$Estimate[plant_div_ci_df$Effect == "(Intercept)"]
intercept_lower <- plant_div_ci_df$CI_lower[plant_div_ci_df$Effect == "(Intercept)"]
intercept_upper <- plant_div_ci_df$CI_upper[plant_div_ci_df$Effect == "(Intercept)"]

# Now create a new column with actual predicted values
plant_div_ci_df$Actual_Estimate <- plant_div_ci_df$Estimate
plant_div_ci_df$Actual_CI_lower <- plant_div_ci_df$CI_lower
plant_div_ci_df$Actual_CI_upper <- plant_div_ci_df$CI_upper

# Add the intercept to all *non-intercept* rows
non_intercept_rows <- plant_div_ci_df$Effect != "(Intercept)"

plant_div_ci_df$Actual_Estimate[non_intercept_rows] <- 
  plant_div_ci_df$Estimate[non_intercept_rows] + intercept_estimate

plant_div_ci_df$Actual_CI_lower[non_intercept_rows] <- 
  plant_div_ci_df$CI_lower[non_intercept_rows] + intercept_lower

plant_div_ci_df$Actual_CI_upper[non_intercept_rows] <- 
  plant_div_ci_df$CI_upper[non_intercept_rows] + intercept_upper


# Using base R
colnames(plant_div_ci_df)[colnames(plant_div_ci_df) == "Effect"] <- "Type"

plant_div_ci_df$Type[plant_div_ci_df$Type == "(Intercept)"] <- "Mature"
plant_div_ci_df$Type[plant_div_ci_df$Type == "TypeRegen"] <- "Regen"
plant_div_ci_df$Type[plant_div_ci_df$Type == "TypeUnforested (heath)"] <- "Unforested (heath)"
plant_div_ci_df$Type[plant_div_ci_df$Type == "TypeUnforested (remnant)"] <- "Unforested (remnant)"

summary(diversity_model)
## Linear mixed-effects model fit by REML
##   Data: plant_div 
##         AIC       BIC   logLik
##   -79.72206 -65.89767 45.86103
## 
## Random effects:
##  Formula: ~1 | Site
##         (Intercept)  Residual
## StdDev: 0.003486986 0.1201534
## 
## Fixed effects:  Diversity_Index ~ Type 
##                               Value  Std.Error DF   t-value p-value
## (Intercept)               0.6861663 0.02629680 72 26.093143  0.0000
## TypeRegen                -0.0356075 0.03708015 72 -0.960284  0.3401
## TypeUnforested (heath)    0.0490610 0.04061927 72  1.207826  0.2311
## TypeUnforested (remnant)  0.0139425 0.03708015 72  0.376009  0.7080
##  Correlation: 
##                          (Intr) TypRgn TypU(h)
## TypeRegen                -0.705               
## TypeUnforested (heath)   -0.644  0.456        
## TypeUnforested (remnant) -0.705  0.500  0.456 
## 
## Standardized Within-Group Residuals:
##        Min         Q1        Med         Q3        Max 
## -3.8956063 -0.5227355  0.1062503  0.6021995  1.6875788 
## 
## Number of Observations: 78
## Number of Groups: 3
plant_div_ci_df
##                                          Type    Estimate    CI_lower
## (Intercept)                            Mature  0.68616626  0.63374454
## TypeRegen                               Regen -0.03560748 -0.10952542
## TypeUnforested (heath)     Unforested (heath)  0.04906101 -0.03191202
## TypeUnforested (remnant) Unforested (remnant)  0.01394246 -0.05997548
##                            CI_upper Actual_Estimate Actual_CI_lower
## (Intercept)              0.73858798       0.6861663       0.6337445
## TypeRegen                0.03831045       0.6505588       0.5242191
## TypeUnforested (heath)   0.13003405       0.7352273       0.6018325
## TypeUnforested (remnant) 0.08786039       0.7001087       0.5737691
##                          Actual_CI_upper
## (Intercept)                    0.7385880
## TypeRegen                      0.7768984
## TypeUnforested (heath)         0.8686220
## TypeUnforested (remnant)       0.8264484

Plotting plant iversity by forest treatment

plant_diversity_graph <- ggplot(plant_div, aes(x = Type, y = Diversity_Index, fill = Type)) +
  
  # Violin plot to show distribution of raw values
  geom_violin(alpha = 0.3, color = NA, width = 0.9) +
  
  # Individual data points (jittered)
  geom_jitter(aes(color = Type), size = 2, shape = 21, width = 0.1, height = 0, alpha = 0.6) +
  
  # Model estimates (mean points from plant_div_ci_df)
  geom_point(data = plant_div_ci_df, 
             aes(x = Type, y = Actual_Estimate), 
             shape = 21, fill = "black", color = "black", size = 3, inherit.aes = FALSE) +
  
  # Confidence intervals (error bars)
  geom_errorbar(data = plant_div_ci_df,
                aes(x = Type, ymin = Actual_CI_lower, ymax = Actual_CI_upper), 
                width = 0.2, color = "black", size = 0.5, inherit.aes = FALSE) +

  # Aesthetics and themes
  theme_bw() +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1),
    strip.background = element_rect(fill = "grey", color = "black", size = 1),
    strip.text = element_text(color = "black", size = 12),
    panel.spacing = unit(1, "lines"),
    panel.background = element_rect(fill = "white", color = "black"),
    plot.margin = margin(10, 10, 10, 10),
    legend.position = "none"
  ) +
  labs(
    x = "Forest Treatment",
    y = "Diversity Index (1 - Simpson's D)",
    fill = "Type",
    title = "Plant Diversity"
  ) +
  scale_fill_manual(values = c(
    "Mature" = "#D55E00",
    "Regen" = "#009E73",
    "Unforested (heath)" = "#56B4E9",
    "Unforested (remnant)" = "#CC79A7"
  )) +
  scale_color_manual(values = c(
    "Mature" = "#D55E00",
    "Regen" = "#009E73",
    "Unforested (heath)" = "#56B4E9",
    "Unforested (remnant)" = "#CC79A7"
  ))

Plant Richness:

Extracting model predicted data:

# Step 1: Extract the intercept value
intercept_value <- coef(summary(richness_model))["(Intercept)", "Estimate"]

# Step 2: Get confidence intervals (for all fixed effects)
wald_ci <- confint(richness_model, level = 0.95)
## Waiting for profiling to be done...
# Step 3: Create results dataframe
results_df <- data.frame(
  Effect = rownames(coef(summary(richness_model))),
  Estimate = coef(summary(richness_model))[, "Estimate"],
  Std_Error = coef(summary(richness_model))[, "Std. Error"],
  p_value = coef(summary(richness_model))[, "Pr(>|z|)"],
  CI_Lower = wald_ci[, 1],
  CI_Upper = wald_ci[, 2]
)

# Step 4: Add real (absolute) estimates and CIs — only modify non-intercepts
results_df$Real_Estimate <- results_df$Estimate
results_df$Real_CI_Lower <- results_df$CI_Lower
results_df$Real_CI_Upper <- results_df$CI_Upper

non_intercepts <- results_df$Effect != "(Intercept)"

results_df$Real_Estimate[non_intercepts] <- results_df$Estimate[non_intercepts] + intercept_value
results_df$Real_CI_Lower[non_intercepts] <- results_df$CI_Lower[non_intercepts] + intercept_value
results_df$Real_CI_Upper[non_intercepts] <- results_df$CI_Upper[non_intercepts] + intercept_value

# Step 5: Rename 'Effect' column to 'Type' and recode factor levels
colnames(results_df)[colnames(results_df) == "Effect"] <- "Type"

results_df$Type[results_df$Type == "(Intercept)"] <- "Mature"
results_df$Type[results_df$Type == "TypeRegen"] <- "Regen"
results_df$Type[results_df$Type == "TypeUnforested (heath)"] <- "Unforested (heath)"
results_df$Type[results_df$Type == "TypeUnforested (remnant)"] <- "Unforested (remnant)"

# Step 6: (Optional) Set row names to match the cleaned-up 'Type' values
rownames(results_df) <- results_df$Type



# Step 8: Exponentiate Estimates and Confidence Intervals
results_df$Real_Estimate <- exp(results_df$Real_Estimate)
results_df$Real_CI_Lower <- exp(results_df$Real_CI_Lower)
results_df$Real_CI_Upper <- exp(results_df$Real_CI_Upper)

print(results_df)
##                                      Type    Estimate  Std_Error       p_value
## Mature                             Mature 2.339972625 0.06772851 1.460262e-261
## Regen                               Regen 0.168464522 0.09199522  6.706602e-02
## Unforested (heath)     Unforested (heath) 0.014572207 0.10448093  8.890769e-01
## Unforested (remnant) Unforested (remnant) 0.009132484 0.09556467  9.238673e-01
##                         CI_Lower  CI_Upper Real_Estimate Real_CI_Lower
## Mature                2.20422005 2.4698410      10.38095      9.063180
## Regen                -0.01150943 0.3493513      12.28571     10.262158
## Unforested (heath)   -0.19150620 0.2184142      10.53333      8.571703
## Unforested (remnant) -0.17828169 0.1966002      10.47619      8.685812
##                      Real_CI_Upper
## Mature                    11.82057
## Regen                     14.72172
## Unforested (heath)        12.91497
## Unforested (remnant)      12.63629

Plotting plant richness by forest treatment

plant_richness_graph <- ggplot(plant_rich, aes(x = Type, y = Richness, fill = Type)) +
  
  # Violin plot to show data distribution
  geom_violin(alpha = 0.3, color = NA, width = 0.9) +

  # Jittered individual data points
  geom_jitter(aes(color = Type), size = 2, shape = 21, width = 0.1, height = 0, alpha = 0.6) +

  # Model estimates (points)
  geom_point(data = results_df, aes(x = Type, y = Real_Estimate),
             shape = 21, fill = "black", size = 3, inherit.aes = FALSE) +

  # Confidence intervals from model
  geom_errorbar(data = results_df, 
                aes(x = Type, ymin = Real_CI_Lower, ymax = Real_CI_Upper),
                width = 0.2, color = "black", size = 0.5, inherit.aes = FALSE) +

  theme_bw() +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1),
    strip.background = element_rect(fill = "grey", color = "black", size = 1),
    strip.text = element_text(color = "black", size = 12),
    panel.spacing = unit(1, "lines"),
    panel.background = element_rect(fill = "white", color = "black"),
    plot.margin = margin(10, 10, 10, 10)
  ) +
  labs(
    x = "Forest Treatment",
    y = "Species Richness",
    fill = "Type",
    title = "Plant Richness"
  ) +
  scale_fill_manual(
    name = "Forest Treatment",
                    values = c(
    "Mature" = "#D55E00",
    "Regen" = "#009E73",
    "Unforested (heath)" = "#56B4E9",
    "Unforested (remnant)" = "#CC79A7"
  )) +
  scale_color_manual(name = "Forest Treatment",
                     values = c(
    "Mature" = "#D55E00",
    "Regen" = "#009E73",
    "Unforested (heath)" = "#56B4E9",
    "Unforested (remnant)" = "#CC79A7"
  ))


plant_rich_div_graph <- plant_diversity_graph + plant_richness_graph

plant_rich_div_graph

save_plot(plant_rich_div_graph, "plant_rich_div_graph")
## png 
##   2

Plant Community Composition:

plant_com_dataset <- plant_com_dataset[-c(561, 701), ]


plant_com_wide <- as.data.frame(spread(plant_com_dataset, Species, DOMIN))

After converting the data into wide format, I need to replace all the NA values with 0

plant_com_wide[is.na(plant_com_wide)] <- 0

plant_com_env <- plant_com_wide[c("Plot", "Site.y", "Type")]

plant_com_wide$Plot <- NULL
plant_com_wide$Site.y <- NULL
plant_com_wide$Type <- NULL
plant_com_wide$Family <- NULL
plant_com_wide$Genus <- NULL

plant_com_wide <- as.data.frame(lapply(plant_com_wide, as.integer))

str(plant_com_wide)

NMDS_plant_com <- metaMDS(plant_com_wide)
NMDS_plant_com

NMDS_plant_com$stress


stressplot(NMDS_plant_com)

plant_com_env$Type <- factor(plant_com_env$Type)


scl <- 4
colvec <- c("#D55E00", "#009E73", "#56B4E9", "#CC79A7")
plot(NMDS_plant_com, type = "n")
with(plant_com_env, points(NMDS_plant_com, display = "sites", col = colvec[Type],
                          scaling = scl, pch = 21, bg = colvec[Type]))
with(plant_com_env, legend("topleft", legend = levels(Type), bty = "n",
                          col = colvec, pch = 21, pt.bg = colvec))

NMDS_plant_com.fit <- envfit(NMDS_plant_com ~ Type , data = plant_com_env, perm = 999)
NMDS_plant_com.fit
## 
## ***FACTORS:
## 
## Centroids:
##                            NMDS1   NMDS2
## TypeMature                0.8232  0.1681
## TypeRegen                 0.0537 -0.1053
## TypeUnforested (heath)   -0.6063 -0.0011
## TypeUnforested (remnant) -0.4438 -0.0620
## 
## Goodness of fit:
##          r2 Pr(>r)    
## Type 0.5383  0.001 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 999
plot(NMDS_plant_com.fit, dis = "site")
plot(NMDS_plant_com.fit)

### Plant NMDS

# 1. Calculate ellipses using vegan::ordiellipse() but don't draw them
library(vegan)

ordi_ellipses <- ordiellipse(NMDS_plant_com, plant_com_env$Type, kind = "se", conf = 0.95, draw = "none")

# 2. Extract ellipse coordinates into a data frame
ellipse_df <- do.call(rbind, lapply(names(ordi_ellipses), function(group) {
  ellipse_points <- vegan:::veganCovEllipse(ordi_ellipses[[group]]$cov, 
                                            ordi_ellipses[[group]]$center, 
                                            ordi_ellipses[[group]]$scale)
  data.frame(ellipse_points, Type = group)
}))
colnames(ellipse_df)[1:2] <- c("NMDS1", "NMDS2")


library(ggplot2)

# Your main NMDS points
nmds_df <- data.frame(
  NMDS1 = NMDS_plant_com$points[, 1],
  NMDS2 = NMDS_plant_com$points[, 2],
  Type = plant_com_env$Type
)

anosim_plant_com_r_tr <- anosim(wisconsin(sqrt(plant_com_wide)), plant_com_env$Type, distance = "bray", permutations = 9999)


R_val <- round(anosim_plant_com_r_tr$statistic, 3)
p_val <- ifelse(anosim_plant_com_r_tr$signif < 0.001, "< 0.001", paste0("= ", signif(anosim_plant_com_r_tr$signif, 2)))

anosim_label <- paste0("ANOSIM R = ", R_val, "\n", "p ", p_val)

# Site scores
site_scores <- as.data.frame(NMDS_plant_com$points)

colnames(site_scores)[1:2] <- c("NMDS1", "NMDS2")

site_scores$Type <- plant_com_env$Type  # Add grouping

# Species scores (usually in $species, or use scores())
species_scores <- as.data.frame(scores(NMDS_plant_com, display = "species"))
species_scores$Species <- rownames(species_scores)

plant_nmds_hull_graph <- ggplot(nmds_df, aes(x = NMDS1, y = NMDS2, color = Type, fill = Type)) +
  geom_point(shape = 21, size = 1, stroke = 1) +
  geom_polygon(data = ellipse_df, aes(x = NMDS1, y = NMDS2, fill = Type, group = Type, color = Type),
               alpha = 0.4, linewidth = 0.5) +
  # ⬇️ Add species labels with inherit.aes = FALSE
  geom_text(data = species_scores, 
            mapping = aes(x = NMDS1, y = NMDS2, label = Species), 
            inherit.aes = FALSE, 
            color = "white", size = 2, alpha = 0) +
  scale_color_manual(name = "Forest Treatment",
                     values = colvec) +
  scale_fill_manual(name = "Forest Treatment",
                    values = colvec) +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
    strip.background = element_rect(fill = "grey", color = "black", size = 1),
    strip.text = element_text(color = "black", size = 12),
    panel.spacing = unit(1, "lines"),
    panel.background = element_rect(fill = "white", color = "black"),
    plot.margin = margin(10, 10, 10, 10),
    legend.position = "none"
) +
  annotate("text", x = Inf, y = Inf, label = anosim_label, hjust = 1.1, vjust = 1.5, size = 3.5, family = "mono") +
  labs(title = "Plant Community Separation");

plant_nmds_hull_graph <- plant_nmds_hull_graph + coord_fixed()

plant_nmds_hull_graph

save_plot(plant_nmds_hull_graph, "plant_nmds_hull_graph")
## png 
##   2

Plant NMDS resampling

# Base ANOSIM R value with all species
base_R <- anosim(wisconsin(sqrt(plant_com_wide)), plant_com_env$Type, distance = "bray", permutations = 999)$statistic

# Store results
jackknife_results <- data.frame(Species = character(),
                                 R_value = numeric(),
                                 Delta_R = numeric(),
                                 stringsAsFactors = FALSE)

# Loop over all species (columns in your wide matrix)
species_names <- colnames(plant_com_wide)

for (sp in species_names) {
  # Remove one species
  reduced_matrix <- plant_com_wide[, setdiff(species_names, sp)]
  
  # Run ANOSIM
  anosim_result <- anosim(wisconsin(sqrt(reduced_matrix)), plant_com_env$Type, distance = "bray", permutations = 999)
  
  # Store R and change from base
  R_val <- anosim_result$statistic
  delta_R <- base_R - R_val
  
  jackknife_results <- rbind(jackknife_results, data.frame(Species = sp, R_value = R_val, Delta_R = delta_R))
}

# Sort by greatest impact (largest Delta_R)
jackknife_results <- jackknife_results[order(-jackknife_results$Delta_R), ]

plant_jk_graph <- ggplot(jackknife_results, aes(x = reorder(Species, Delta_R), y = Delta_R)) +
  geom_col(fill = "#0072B2") +
  coord_flip() +
  labs(title = "Plant Species Impact on NMDS Clustering",
       x = "Species Removed",
       y = "Decrease in ANOSIM R") +
  theme_minimal()

plant_jk_graph

save_plot(plant_jk_graph, "plant_jk")
## png 
##   2

Plant NMD (No Trees)

# Define species to exclude
tree_species <- c("Silver birch", "Rowan", "Downy birch")

# Filter matrix to exclude these species
plant_com_no_trees <- plant_com_wide[, !colnames(plant_com_wide) %in% tree_species]
# Run NMDS on the modified dataset
NMDS_plant_com_no_trees <- metaMDS(wisconsin(sqrt(plant_com_no_trees)), distance = "bray", k = 2, trymax = 100)
## Run 0 stress 0.1901061 
## Run 1 stress 0.21296 
## Run 2 stress 0.1883802 
## ... New best solution
## ... Procrustes: rmse 0.03239242  max resid 0.1414564 
## Run 3 stress 0.1890093 
## Run 4 stress 0.194464 
## Run 5 stress 0.1930042 
## Run 6 stress 0.2142864 
## Run 7 stress 0.2104666 
## Run 8 stress 0.214057 
## Run 9 stress 0.2076797 
## Run 10 stress 0.1977019 
## Run 11 stress 0.2114374 
## Run 12 stress 0.2116637 
## Run 13 stress 0.2044696 
## Run 14 stress 0.196647 
## Run 15 stress 0.1983483 
## Run 16 stress 0.188837 
## ... Procrustes: rmse 0.01125132  max resid 0.08060584 
## Run 17 stress 0.2083896 
## Run 18 stress 0.2083942 
## Run 19 stress 0.2040663 
## Run 20 stress 0.2083267 
## Run 21 stress 0.1881442 
## ... New best solution
## ... Procrustes: rmse 0.02078665  max resid 0.1242391 
## Run 22 stress 0.2101414 
## Run 23 stress 0.2157395 
## Run 24 stress 0.2065798 
## Run 25 stress 0.2136877 
## Run 26 stress 0.1972916 
## Run 27 stress 0.197115 
## Run 28 stress 0.1943633 
## Run 29 stress 0.2093802 
## Run 30 stress 0.1890592 
## Run 31 stress 0.1986611 
## Run 32 stress 0.2160472 
## Run 33 stress 0.2086404 
## Run 34 stress 0.2102971 
## Run 35 stress 0.1887508 
## Run 36 stress 0.2195215 
## Run 37 stress 0.1903416 
## Run 38 stress 0.188585 
## ... Procrustes: rmse 0.01487024  max resid 0.1208851 
## Run 39 stress 0.1971506 
## Run 40 stress 0.2079329 
## Run 41 stress 0.1928808 
## Run 42 stress 0.2134388 
## Run 43 stress 0.219096 
## Run 44 stress 0.2132646 
## Run 45 stress 0.2222949 
## Run 46 stress 0.198661 
## Run 47 stress 0.2145521 
## Run 48 stress 0.215177 
## Run 49 stress 0.214961 
## Run 50 stress 0.1881211 
## ... New best solution
## ... Procrustes: rmse 0.001937691  max resid 0.01359463 
## Run 51 stress 0.1900061 
## Run 52 stress 0.2081687 
## Run 53 stress 0.2173113 
## Run 54 stress 0.1974102 
## Run 55 stress 0.2001864 
## Run 56 stress 0.2044841 
## Run 57 stress 0.1961687 
## Run 58 stress 0.2200975 
## Run 59 stress 0.2147783 
## Run 60 stress 0.2109233 
## Run 61 stress 0.2042485 
## Run 62 stress 0.2134556 
## Run 63 stress 0.1941284 
## Run 64 stress 0.1925906 
## Run 65 stress 0.2060116 
## Run 66 stress 0.1886733 
## Run 67 stress 0.189057 
## Run 68 stress 0.2013378 
## Run 69 stress 0.2017538 
## Run 70 stress 0.1881142 
## ... New best solution
## ... Procrustes: rmse 0.002560617  max resid 0.01920495 
## Run 71 stress 0.1885752 
## ... Procrustes: rmse 0.01453685  max resid 0.1212143 
## Run 72 stress 0.2116655 
## Run 73 stress 0.2082069 
## Run 74 stress 0.2039048 
## Run 75 stress 0.1888904 
## Run 76 stress 0.18963 
## Run 77 stress 0.2034862 
## Run 78 stress 0.2066994 
## Run 79 stress 0.1954976 
## Run 80 stress 0.2052672 
## Run 81 stress 0.222485 
## Run 82 stress 0.189617 
## Run 83 stress 0.2036731 
## Run 84 stress 0.2016689 
## Run 85 stress 0.1884583 
## ... Procrustes: rmse 0.02027165  max resid 0.1223164 
## Run 86 stress 0.2166119 
## Run 87 stress 0.2075466 
## Run 88 stress 0.2079606 
## Run 89 stress 0.1886734 
## Run 90 stress 0.2018272 
## Run 91 stress 0.1966538 
## Run 92 stress 0.2178612 
## Run 93 stress 0.1888607 
## Run 94 stress 0.1977504 
## Run 95 stress 0.2016763 
## Run 96 stress 0.1881814 
## ... Procrustes: rmse 0.0319324  max resid 0.1595439 
## Run 97 stress 0.2066725 
## Run 98 stress 0.2061703 
## Run 99 stress 0.2073608 
## Run 100 stress 0.1932321 
## *** Best solution was not repeated -- monoMDS stopping criteria:
##     97: stress ratio > sratmax
##      3: scale factor of the gradient < sfgrmin
anosim_no_trees <- anosim(wisconsin(sqrt(plant_com_no_trees)), plant_com_env$Type, distance = "bray", permutations = 9999)

# Create label for the plot
R_val <- round(anosim_no_trees$statistic, 3)
p_val <- ifelse(anosim_no_trees$signif < 0.001, "< 0.001", paste0("= ", signif(anosim_no_trees$signif, 2)))
anosim_label <- paste0("ANOSIM R = ", R_val, "\n", "p ", p_val)
ordi_ellipses_no_trees <- ordiellipse(NMDS_plant_com_no_trees, plant_com_env$Type, kind = "se", conf = 0.95, draw = "none")

ellipse_df_no_trees <- do.call(rbind, lapply(names(ordi_ellipses_no_trees), function(group) {
  ellipse_points <- vegan:::veganCovEllipse(ordi_ellipses_no_trees[[group]]$cov, 
                                            ordi_ellipses_no_trees[[group]]$center, 
                                            ordi_ellipses_no_trees[[group]]$scale)
  data.frame(ellipse_points, Type = group)
}))
colnames(ellipse_df_no_trees)[1:2] <- c("NMDS1", "NMDS2")
# Site scores
nmds_df_no_trees <- data.frame(
  NMDS1 = NMDS_plant_com_no_trees$points[, 1],
  NMDS2 = NMDS_plant_com_no_trees$points[, 2],
  Type = plant_com_env$Type
)

# Species scores for labels (optional, may require handling NAs)
species_scores <- data.frame(scores(NMDS_plant_com_no_trees, display = "species"))
species_scores$Species <- rownames(species_scores)
colnames(species_scores)[1:2] <- c("NMDS1", "NMDS2")
plant_nmds_no_trees_graph <- ggplot(nmds_df_no_trees, aes(x = NMDS1, y = NMDS2, color = Type, fill = Type)) +
  geom_point(shape = 21, size = 1, stroke = 1) +
  geom_polygon(data = ellipse_df_no_trees, aes(x = NMDS1, y = NMDS2, fill = Type, group = Type, color = Type),
               alpha = 0.4, linewidth = 0.5) +
  geom_text(data = species_scores, 
            aes(x = NMDS1, y = NMDS2, label = Species), 
            inherit.aes = FALSE, 
            color = "white", size = 2, alpha = 0) +
  scale_color_manual(name = "Forest Treatment",
                     values = colvec) +
  scale_fill_manual(name = "Forest Treatment",
                    values = colvec) +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
    strip.background = element_rect(fill = "grey", color = "black", size = 1),
    strip.text = element_text(color = "black", size = 12),
    panel.spacing = unit(1, "lines"),
    panel.background = element_rect(fill = "white", color = "black"),
    plot.margin = margin(10, 10, 10, 10),) +
  annotate("text", x = Inf, y = Inf, label = anosim_label, hjust = 1.1, vjust = 1.5, size = 3.5, family = "mono") +
  labs(title = "Plant Community Separation (Without Tree Species)")

plant_nmds_no_trees_graph <- plant_nmds_no_trees_graph + coord_fixed()

plant_nmds_no_trees_graph

save_plot(plant_nmds_no_trees_graph, "plant_nmds_no_trees_graph")
## png 
##   2

Plant NMDS resampling (No Trees)

# Base ANOSIM R value with no tree species
base_R_no_trees <- anosim(wisconsin(sqrt(plant_com_no_trees)), plant_com_env$Type, distance = "bray", permutations = 999)$statistic

# Store results
jackknife_results_no_trees <- data.frame(Species = character(),
                                         R_value = numeric(),
                                         Delta_R = numeric(),
                                         stringsAsFactors = FALSE)

# Loop over all species (columns) in the no-trees matrix
species_names_no_trees <- colnames(plant_com_no_trees)

for (sp in species_names_no_trees) {
  # Remove one species
  reduced_matrix <- plant_com_no_trees[, setdiff(species_names_no_trees, sp)]
  
  # Run ANOSIM
  anosim_result <- anosim(wisconsin(sqrt(reduced_matrix)), plant_com_env$Type, distance = "bray", permutations = 999)
  
  # Store R and change from base
  R_val <- anosim_result$statistic
  delta_R <- base_R_no_trees - R_val
  
  jackknife_results_no_trees <- rbind(jackknife_results_no_trees,
                                      data.frame(Species = sp, R_value = R_val, Delta_R = delta_R))
}

# Sort by greatest impact (largest Delta_R)
jackknife_results_no_trees <- jackknife_results_no_trees[order(-jackknife_results_no_trees$Delta_R), ]

# Plot the results
plant_jk_no_trees_graph <- ggplot(jackknife_results_no_trees, aes(x = reorder(Species, Delta_R), y = Delta_R)) +
  geom_col(fill = "#0072B2") +
  coord_flip() +
  labs(title = "Plant Species Impact on NMDS Clustering (No Trees)",
       x = "Species Removed",
       y = "Decrease in ANOSIM R") +
  theme_minimal()

jackknife_results_no_trees
##                      Species   R_value       Delta_R
## 44              Silver.birch 0.3659864  4.875403e-02
## 24      Hylocomium.splendens 0.3755024  3.923802e-02
## 3                  Blaeberry 0.3780339  3.670650e-02
## 17               Downy.birch 0.3845283  3.021212e-02
## 2                   Bent.sp. 0.3984823  1.625814e-02
## 10                  Cowberry 0.4054596  9.280855e-03
## 12      Cross.leaved.heather 0.4090427  5.697729e-03
## 20            Heath.bedstraw 0.4107332  4.007247e-03
## 5                 Bog.myrtle 0.4108004  3.940059e-03
## 33        Pleurozia.purpurea 0.4112755  3.464949e-03
## 19                 Hard.fern 0.4120703  2.670098e-03
## 4               Bog.asphodel 0.4124495  2.290969e-03
## 15      Devil.s.bit.scabious 0.4127284  2.012022e-03
## 39    Rhytidiadelphus.loreus 0.4133907  1.349746e-03
## 36           Polytrichum.sp. 0.4136990  1.041404e-03
## 47                Star.sedge 0.4137026  1.037805e-03
## 42                Scots.pine 0.4137062  1.034206e-03
## 1               Bell.heather 0.4137752  9.652185e-04
## 13                 Crowberry 0.4137758  9.646187e-04
## 32                  Moss.sp. 0.4139216  8.188461e-04
## 9           Cotton.grass.sp. 0.4139330  8.074482e-04
## 21            Heath.milkwort 0.4140314  7.090667e-04
## 16    Devil.s.bit.scabious.1 0.4140961  6.442789e-04
## 41       Round.leaved.sundew 0.4143157  4.247202e-04
## 38           Racomitrium.sp. 0.4144393  3.011434e-04
## 26              Ling.heather 0.4146858  5.458974e-05
## 52        Unknown.moss.sp..3 0.4147848 -4.439165e-05
## 7                  Club.moss 0.4148028 -6.238827e-05
## 28           Marsh.lousewort 0.4149264 -1.859650e-04
## 54              Wood.anemone 0.4151759 -4.355181e-04
## 11                   Cowslip 0.4152011 -4.607134e-04
## 56              Woodrush.sp. 0.4154225 -6.820718e-04
## 27             Liverwort.sp. 0.4157674 -1.027007e-03
## 25        Lemon.scented.fern 0.4158274 -1.086996e-03
## 8           Common.cow.wheat 0.4158874 -1.146984e-03
## 37 Ptilium.crista.castrensis 0.4162186 -1.478122e-03
## 51        Unknown.moss.sp..2 0.4162515 -1.511116e-03
## 18                 Eyebright 0.4163085 -1.568105e-03
## 40          Ribwort.plantain 0.4165155 -1.775066e-03
## 23             Horsetail.sp. 0.4165377 -1.797262e-03
## 49          Tufted.hairgrass 0.4165377 -1.797262e-03
## 22           Heath.speedwell 0.4167656 -2.025219e-03
## 29          Meadow.buttercup 0.4168730 -2.132599e-03
## 35                   Poa.sp. 0.4169684 -2.227981e-03
## 55               Wood.sorrel 0.4170224 -2.281971e-03
## 30          Meadow.grass.sp. 0.4170872 -2.346759e-03
## 6                    Bracken 0.4172384 -2.497930e-03
## 14                 Deergrass 0.4175659 -2.825469e-03
## 50          Unknown.moss.sp. 0.4178514 -3.111015e-03
## 45                 Soft.rush 0.4180110 -3.270585e-03
## 43                 Sedge.sp. 0.4181508 -3.410359e-03
## 53                Violet.sp. 0.4189186 -4.178214e-03
## 34      Pleurozium.schreberi 0.4198161 -5.075646e-03
## 31          Molinia.caerulea 0.4229829 -8.242450e-03
## 48                 Tormentil 0.4470726 -3.233212e-02
## 46              Sphagnum.sp. 0.4488872 -3.414678e-02
# Display the graph
plant_jk_no_trees_graph

save_plot(plant_jk_no_trees_graph, "plant_jk_no_trees")
## png 
##   2
palette.colors()
## [1] "#000000" "#E69F00" "#56B4E9" "#009E73" "#F0E442" "#0072B2" "#D55E00"
## [8] "#CC79A7" "#999999"

Soil Invertebrate Diversity:

Extracting model predicted data

predicted_soil_inv_diversity <- predict(soil_invert_diversity_model, type = "response", level = 1)

predicted_soil_div_with_ci <- intervals(soil_invert_diversity_model, level = 0.95)

# Extract just the fixed effects into a data frame
soil_div_ci_df <- as.data.frame(predicted_soil_div_with_ci$fixed)


# Optionally, name the columns
colnames(soil_div_ci_df) <- c("CI_lower", "Estimate", "CI_upper")

# Add a column for the effect names (e.g., Intercept, Carbon, etc.)
soil_div_ci_df$Effect <- rownames(predicted_soil_div_with_ci$fixed)

# Reorder for readability
soil_div_ci_df <- soil_div_ci_df[, c("Effect", "Estimate", "CI_lower", "CI_upper")]

# View it
print(soil_div_ci_df)
##                                            Effect    Estimate    CI_lower
## (Intercept)                           (Intercept)  0.36597967  0.26761870
## TypeRegen                               TypeRegen  0.09497504 -0.04674363
## TypeUnforested (heath)     TypeUnforested (heath) -0.06544998 -0.21779287
## TypeUnforested (remnant) TypeUnforested (remnant) -0.05195028 -0.19366895
##                            CI_upper
## (Intercept)              0.46434064
## TypeRegen                0.23669371
## TypeUnforested (heath)   0.08689291
## TypeUnforested (remnant) 0.08976839
# Extract the intercept value
soil_intercept_estimate <- soil_div_ci_df$Estimate[soil_div_ci_df$Effect == "(Intercept)"]
soil_intercept_lower <- soil_div_ci_df$CI_lower[soil_div_ci_df$Effect == "(Intercept)"]
soil_intercept_upper <- soil_div_ci_df$CI_upper[soil_div_ci_df$Effect == "(Intercept)"]

# Now create a new column with actual predicted values
soil_div_ci_df$Actual_Estimate <- soil_div_ci_df$Estimate
soil_div_ci_df$Actual_CI_lower <- soil_div_ci_df$CI_lower
soil_div_ci_df$Actual_CI_upper <- soil_div_ci_df$CI_upper

# Add the intercept to all *non-intercept* rows
soil_non_intercept_rows <- soil_div_ci_df$Effect != "(Intercept)"

soil_div_ci_df$Actual_Estimate[soil_non_intercept_rows] <- 
  soil_div_ci_df$Estimate[soil_non_intercept_rows] + soil_intercept_estimate

soil_div_ci_df$Actual_CI_lower[soil_non_intercept_rows] <- 
  soil_div_ci_df$CI_lower[soil_non_intercept_rows] + soil_intercept_lower

soil_div_ci_df$Actual_CI_upper[soil_non_intercept_rows] <- 
  soil_div_ci_df$CI_upper[soil_non_intercept_rows] + soil_intercept_upper


# Using base R
colnames(soil_div_ci_df)[colnames(soil_div_ci_df) == "Effect"] <- "Type"

soil_div_ci_df$Type[soil_div_ci_df$Type == "(Intercept)"] <- "Mature"
soil_div_ci_df$Type[soil_div_ci_df$Type == "TypeRegen"] <- "Regen"
soil_div_ci_df$Type[soil_div_ci_df$Type == "TypeUnforested (heath)"] <- "Unforested (heath)"
soil_div_ci_df$Type[soil_div_ci_df$Type == "TypeUnforested (remnant)"] <- "Unforested (remnant)"


soil_div_ci_df
##                                          Type    Estimate    CI_lower
## (Intercept)                            Mature  0.36597967  0.26761870
## TypeRegen                               Regen  0.09497504 -0.04674363
## TypeUnforested (heath)     Unforested (heath) -0.06544998 -0.21779287
## TypeUnforested (remnant) Unforested (remnant) -0.05195028 -0.19366895
##                            CI_upper Actual_Estimate Actual_CI_lower
## (Intercept)              0.46434064       0.3659797      0.26761870
## TypeRegen                0.23669371       0.4609547      0.22087506
## TypeUnforested (heath)   0.08689291       0.3005297      0.04982583
## TypeUnforested (remnant) 0.08976839       0.3140294      0.07394974
##                          Actual_CI_upper
## (Intercept)                    0.4643406
## TypeRegen                      0.7010343
## TypeUnforested (heath)         0.5512335
## TypeUnforested (remnant)       0.5541090

Plotting soil invertebrate diversity

soil_diversity_graph <- ggplot(soil_invert_div, aes(x = Type, y = Diversity_Index, fill = Type)) +
  
  # Violin plot to show distribution of raw diversity index values
  geom_violin(alpha = 0.3, color = NA, width = 0.9) +

  # Individual data points (jittered)
  geom_jitter(aes(color = Type), size = 2, shape = 21, width = 0.1, height = 0, alpha = 0.6) +

  # Model-estimated means (points)
  geom_point(data = soil_div_ci_df, 
             aes(x = Type, y = Actual_Estimate), 
             shape = 21, fill = "black", color = "black", size = 3, inherit.aes = FALSE) +

  # Confidence intervals (error bars)
  geom_errorbar(data = soil_div_ci_df, 
                aes(x = Type, ymin = Actual_CI_lower, ymax = Actual_CI_upper), 
                width = 0.2, color = "black", size = 0.5, inherit.aes = FALSE) +

  # Theme and styling
  theme_bw() +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1),
    strip.background = element_rect(fill = "grey", color = "black", size = 1),
    strip.text = element_text(color = "black", size = 12),
    panel.spacing = unit(1, "lines"),
    panel.background = element_rect(fill = "white", color = "black"),
    plot.margin = margin(10, 10, 10, 10),
    legend.position = "none"
  ) +
  labs(
    x = "Forest Treatment",
    y = "Diversity Index (1 - Simpson's D)",
    fill = "Type",
    title = "Soil Invertebrate Diversity"
  ) +
  scale_fill_manual(
    name = "Forest Treatment",
    values = c(
    "Mature" = "#D55E00",
    "Regen" = "#009E73",
    "Unforested (heath)" = "#56B4E9",
    "Unforested (remnant)" = "#CC79A7"
  )) +
  scale_color_manual(
    name = "Forest Treatment",
    values = c(
    "Mature" = "#D55E00",
    "Regen" = "#009E73",
    "Unforested (heath)" = "#56B4E9",
    "Unforested (remnant)" = "#CC79A7"
  ))

Soil Invertebrate Richness:

Extracting soil invertebrate richness data:

# Step 1: Extract the intercept value
soil_intercept_value <- coef(summary(soil_invert_richness_model))["(Intercept)", "Estimate"]

soil_wald_ci <- confint(soil_invert_richness_model, level = 0.95)
## Waiting for profiling to be done...
# Drop .sig01 row
soil_wald_ci <- soil_wald_ci[!rownames(soil_wald_ci) %in% ".sig01", ]

# Then build the dataframe
soil_results_df <- data.frame(
  Effect = rownames(coef(summary(soil_invert_richness_model))),
  Estimate = coef(summary(soil_invert_richness_model))[, "Estimate"],
  Std_Error = coef(summary(soil_invert_richness_model))[, "Std. Error"],
  p_value = coef(summary(soil_invert_richness_model))[, "Pr(>|z|)"],
  CI_Lower = soil_wald_ci[, 1],
  CI_Upper = soil_wald_ci[, 2]
)


# Step 4: Add real (absolute) estimates and CIs — only modify non-intercepts
soil_results_df$Real_Estimate <- soil_results_df$Estimate
soil_results_df$Real_CI_Lower <- soil_results_df$CI_Lower
soil_results_df$Real_CI_Upper <- soil_results_df$CI_Upper

soil_non_intercepts <- soil_results_df$Effect != "(Intercept)"

soil_results_df$Real_Estimate[soil_non_intercepts] <- soil_results_df$Estimate[soil_non_intercepts] + soil_intercept_value
soil_results_df$Real_CI_Lower[soil_non_intercepts] <- soil_results_df$CI_Lower[soil_non_intercepts] + soil_intercept_value
soil_results_df$Real_CI_Upper[soil_non_intercepts] <- soil_results_df$CI_Upper[soil_non_intercepts] + soil_intercept_value

# Step 5: Rename 'Effect' column to 'Type' and recode factor levels
colnames(soil_results_df)[colnames(soil_results_df) == "Effect"] <- "Type"

soil_results_df$Type[soil_results_df$Type == "(Intercept)"] <- "Mature"
soil_results_df$Type[soil_results_df$Type == "TypeRegen"] <- "Regen"
soil_results_df$Type[soil_results_df$Type == "TypeUnforested (heath)"] <- "Unforested (heath)"
soil_results_df$Type[soil_results_df$Type == "TypeUnforested (remnant)"] <- "Unforested (remnant)"

# Step 6: (Optional) Set row names to match the cleaned-up 'Type' values
rownames(soil_results_df) <- soil_results_df$Type



# Step 8: Exponentiate Estimates and Confidence Intervals
soil_results_df$Real_Estimate <- exp(soil_results_df$Real_Estimate)
soil_results_df$Real_CI_Lower <- exp(soil_results_df$Real_CI_Lower)
soil_results_df$Real_CI_Upper <- exp(soil_results_df$Real_CI_Upper)

print(soil_results_df)
##                                      Type   Estimate Std_Error      p_value
## Mature                             Mature  1.3312346 0.1373606 3.275698e-22
## Regen                               Regen -0.2852660 0.2142310 1.829978e-01
## Unforested (heath)     Unforested (heath) -0.3379828 0.2364418 1.528736e-01
## Unforested (remnant) Unforested (remnant) -0.3996764 0.2217452 7.148053e-02
##                        CI_Lower   CI_Upper Real_Estimate Real_CI_Lower
## Mature                1.0493695 1.58889621      3.785714      2.855850
## Regen                -0.7119726 0.13093157      2.846154      1.857557
## Unforested (heath)   -0.8152634 0.11588244      2.700000      1.675265
## Unforested (remnant) -0.8434047 0.02932028      2.538462      1.628778
##                      Real_CI_Upper
## Mature                    4.898339
## Regen                     4.315297
## Unforested (heath)        4.250842
## Unforested (remnant)      3.898356

Plotting soil invertebrate richness

soil_invert_richness_graph <- ggplot(soil_invert_rich, aes(x = Type, y = Richness, fill = Type)) +

  # Violin plot for distribution of raw richness values
  geom_violin(alpha = 0.3, color = NA, width = 0.9) +

  # Jittered individual data points
  geom_jitter(aes(color = Type), size = 2, shape = 21, width = 0.1, height = 0, alpha = 0.6) +

  # Model-based means (points)
  geom_point(data = soil_results_df,
             aes(x = Type, y = Real_Estimate),
             shape = 21, fill = "black", color = "black", size = 3, inherit.aes = FALSE) +

  # Model-based confidence intervals (error bars)
  geom_errorbar(data = soil_results_df,
                aes(x = Type, ymin = Real_CI_Lower, ymax = Real_CI_Upper),
                width = 0.2, color = "black", size = 0.5, inherit.aes = FALSE) +

  # Clean minimal theme
  theme_bw() +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1),
    strip.background = element_rect(fill = "grey", color = "black", size = 1),
    strip.text = element_text(color = "black", size = 12),
    panel.spacing = unit(1, "lines"),
    panel.background = element_rect(fill = "white", color = "black"),
    plot.margin = margin(10, 10, 10, 10)
  ) +
  labs(
    x = "Forest Treatment",
    y = "Order level Richness",
    fill = "Type",
    title = "Soil invertebrate Richness"
  ) +
  scale_fill_manual(
    name = "Forest Treatment",
    values = c(
    "Mature" = "#D55E00",
    "Regen" = "#009E73",
    "Unforested (heath)" = "#56B4E9",
    "Unforested (remnant)" = "#CC79A7"
  )) +
  scale_color_manual(
    name = "Forest Treatment",
    values = c(
    "Mature" = "#D55E00",
    "Regen" = "#009E73",
    "Unforested (heath)" = "#56B4E9",
    "Unforested (remnant)" = "#CC79A7"
  ))

soil_div_rich_graph <- soil_diversity_graph + soil_invert_richness_graph

soil_div_rich_graph

save_plot(soil_div_rich_graph, "soil_rich_div_graph")
## png 
##   2

Soil Invertebrate Community Composition:

#soil_invert_dataset <- plant_com_dataset[-c(561, 701), ]

soil_invert_dataset <- soil_invert_dataset[c("Order", "Type", "Abundance", "Site.y", "Plot")]


soil_invert_dataset <- aggregate(cbind(Abundance) ~ Order + Plot, data = soil_invert_dataset, FUN = sum)

soil_invert_dataset <- merge(soil_invert_dataset, overview_dataset, by = "Plot")

soil_invert_dataset <- soil_invert_dataset[c("Order", "Type", "Abundance", "Site", "Plot")]

soil_invert_wide <- as.data.frame(spread(soil_invert_dataset, Order, Abundance))

After converting the data into wide format, I need to replace all the NA values with 0

soil_invert_wide[is.na(soil_invert_wide)] <- 0

soil_invert_env <- soil_invert_wide[c("Plot", "Type", "Site")]

soil_invert_wide$Plot <- NULL
soil_invert_wide$Site <- NULL
soil_invert_wide$Type <- NULL
soil_invert_wide$Family <- NULL
soil_invert_wide$Genus <- NULL

soil_invert_wide <- as.data.frame(lapply(soil_invert_wide, as.integer))

str(soil_invert_wide)
## 'data.frame':    50 obs. of  11 variables:
##  $ Acari           : int  11 14 217 75 304 10 55 70 109 27 ...
##  $ Araneae         : int  0 0 0 0 0 0 0 0 0 1 ...
##  $ Coleoptera      : int  0 0 2 0 0 0 0 0 0 0 ...
##  $ Diptera         : int  4 1 2 1 4 5 4 2 1 8 ...
##  $ Entomobryomorpha: int  0 0 0 0 0 1 0 0 0 0 ...
##  $ Geophilomorpha  : int  0 0 0 0 0 0 0 1 0 0 ...
##  $ Hemiptera       : int  0 0 0 0 10 0 0 0 0 0 ...
##  $ Hymenoptera     : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Isopoda         : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Podomorpha      : int  5 1 70 9 46 0 7 61 9 3 ...
##  $ Thysanoptera    : int  0 0 0 0 0 0 0 0 0 0 ...
NMDS_soil_invert <- metaMDS(soil_invert_wide)
## Square root transformation
## Wisconsin double standardization
## Run 0 stress 0.1761955 
## Run 1 stress 0.1712266 
## ... New best solution
## ... Procrustes: rmse 0.04440726  max resid 0.2888872 
## Run 2 stress 0.2220741 
## Run 3 stress 0.1926232 
## Run 4 stress 0.2231124 
## Run 5 stress 0.1725766 
## Run 6 stress 0.1836321 
## Run 7 stress 0.1843298 
## Run 8 stress 0.176317 
## Run 9 stress 0.1730456 
## Run 10 stress 0.197426 
## Run 11 stress 0.1884122 
## Run 12 stress 0.2083265 
## Run 13 stress 0.1992799 
## Run 14 stress 0.1906248 
## Run 15 stress 0.2098807 
## Run 16 stress 0.1727093 
## Run 17 stress 0.1744323 
## Run 18 stress 0.1773268 
## Run 19 stress 0.1758481 
## Run 20 stress 0.1776654 
## *** Best solution was not repeated -- monoMDS stopping criteria:
##      1: no. of iterations >= maxit
##     19: stress ratio > sratmax
NMDS_soil_invert
## 
## Call:
## metaMDS(comm = soil_invert_wide) 
## 
## global Multidimensional Scaling using monoMDS
## 
## Data:     wisconsin(sqrt(soil_invert_wide)) 
## Distance: bray 
## 
## Dimensions: 2 
## Stress:     0.1712266 
## Stress type 1, weak ties
## Best solution was not repeated after 20 tries
## The best solution was from try 1 (random start)
## Scaling: centring, PC rotation, halfchange scaling 
## Species: expanded scores based on 'wisconsin(sqrt(soil_invert_wide))'
stressplot(NMDS_soil_invert)

soil_invert_env
##    Plot                 Type         Site
## 1    63               Mature Glen Cluanie
## 2    66               Mature Glen Cluanie
## 3    68               Mature Glen Cluanie
## 4    69               Mature Glen Cluanie
## 5    74               Mature Glen Cluanie
## 6    81               Mature Glen Cluanie
## 7    82               Mature Glen Cluanie
## 8     3               Mature  Glen Mallie
## 9     4               Mature  Glen Mallie
## 10   10               Mature  Glen Mallie
## 11   12               Mature  Glen Mallie
## 12   14               Mature  Glen Mallie
## 13   24               Mature  Glen Mallie
## 14   25               Mature  Glen Mallie
## 15   62                Regen Glen Cluanie
## 16   64                Regen Glen Cluanie
## 17   65                Regen Glen Cluanie
## 18   67                Regen Glen Cluanie
## 19   71                Regen Glen Cluanie
## 20   80                Regen Glen Cluanie
## 21    6                Regen  Glen Mallie
## 22    7                Regen  Glen Mallie
## 23    8                Regen  Glen Mallie
## 24   11                Regen  Glen Mallie
## 25   13                Regen  Glen Mallie
## 26   23                Regen  Glen Mallie
## 27   26                Regen  Glen Mallie
## 28   57   Unforested (heath) Glen Cluanie
## 29   58   Unforested (heath) Glen Cluanie
## 30   59   Unforested (heath) Glen Cluanie
## 31   60   Unforested (heath) Glen Cluanie
## 32   77   Unforested (heath) Glen Cluanie
## 33   16   Unforested (heath)  Glen Mallie
## 34   17   Unforested (heath)  Glen Mallie
## 35   18   Unforested (heath)  Glen Mallie
## 36   19   Unforested (heath)  Glen Mallie
## 37   20   Unforested (heath)  Glen Mallie
## 38   61 Unforested (remnant) Glen Cluanie
## 39   70 Unforested (remnant) Glen Cluanie
## 40   72 Unforested (remnant) Glen Cluanie
## 41   75 Unforested (remnant) Glen Cluanie
## 42   76 Unforested (remnant) Glen Cluanie
## 43   78 Unforested (remnant) Glen Cluanie
## 44    1 Unforested (remnant)  Glen Mallie
## 45    2 Unforested (remnant)  Glen Mallie
## 46    5 Unforested (remnant)  Glen Mallie
## 47    9 Unforested (remnant)  Glen Mallie
## 48   15 Unforested (remnant)  Glen Mallie
## 49   22 Unforested (remnant)  Glen Mallie
## 50   27 Unforested (remnant)  Glen Mallie
soil_invert_env$Type <- factor(soil_invert_env$Type)



scl <- 4
colvec <- c("#D55E00", "#009E73", "#56B4E9", "#CC79A7")
plot(NMDS_soil_invert, type = "n")
with(soil_invert_env, points(NMDS_soil_invert, display = "sites", col = colvec[Type],
                          scaling = scl, pch = 21, bg = colvec[Type]))
## Warning in plot.xy(xy.coords(x, y), type = type, ...): "scaling" is not a
## graphical parameter
with(soil_invert_env, legend("topleft", legend = levels(Type), bty = "n",
                          col = colvec, pch = 21, pt.bg = colvec))

NMDS_soil_invert.fit <- envfit(NMDS_soil_invert ~ Type , data = soil_invert_env, perm = 999)
NMDS_soil_invert.fit
## 
## ***FACTORS:
## 
## Centroids:
##                            NMDS1   NMDS2
## TypeMature               -0.0571  0.2111
## TypeRegen                -0.0450  0.1406
## TypeUnforested (heath)   -0.0242 -0.1613
## TypeUnforested (remnant)  0.1251 -0.2439
## 
## Goodness of fit:
##          r2 Pr(>r)
## Type 0.0823  0.224
## Permutation: free
## Number of permutations: 999
plot(NMDS_soil_invert.fit, dis = "site")
## Warning in text.default(x$factors$centroids[, choices, drop = FALSE], labs$f, :
## "dis" is not a graphical parameter
plot(NMDS_soil_invert.fit)

### Soil invertebrate NMDS

# 1. Calculate ellipses using vegan::ordiellipse() but don't draw them
library(vegan)

soil_ordi_ellipses <- ordiellipse(NMDS_soil_invert, soil_invert_env$Type, kind = "se", conf = 0.95, draw = "none")

# 2. Extract ellipse coordinates into a data frame
soil_ellipse_df <- do.call(rbind, lapply(names(soil_ordi_ellipses), function(group) {
  soil_ellipse_points <- vegan:::veganCovEllipse(soil_ordi_ellipses[[group]]$cov, 
                                            soil_ordi_ellipses[[group]]$center, 
                                            soil_ordi_ellipses[[group]]$scale)
  data.frame(soil_ellipse_points, Type = group)
}))
colnames(soil_ellipse_df)[1:2] <- c("NMDS1", "NMDS2")


library(ggplot2)

# Your main NMDS points
soil_nmds_df <- data.frame(
  NMDS1 = NMDS_soil_invert$points[, 1],
  NMDS2 = NMDS_soil_invert$points[, 2],
  Type = soil_invert_env$Type
)

soil_anosim_plant_com_r_tr <- anosim(wisconsin(sqrt(soil_invert_wide)), soil_invert_env$Type, distance = "bray", permutations = 9999)


soil_R_val <- round(soil_anosim_plant_com_r_tr$statistic, 3)
soil_p_val <- ifelse(soil_anosim_plant_com_r_tr$signif < 0.001, "< 0.001", paste0("= ", signif(soil_anosim_plant_com_r_tr$signif, 2)))

soil_anosim_label <- paste0("ANOSIM R = ", soil_R_val, "\n", "p ", soil_p_val)

# Site scores
soil_site_scores <- as.data.frame(NMDS_soil_invert$points)

colnames(soil_site_scores)[1:2] <- c("NMDS1", "NMDS2")

soil_site_scores$Type <- soil_invert_env$Type  # Add grouping

# Species scores (usually in $species, or use scores())
soil_species_scores <- as.data.frame(scores(NMDS_soil_invert, display = "species"))
soil_species_scores$Species <- rownames(soil_species_scores)

soil_nmds_hull_graph <- ggplot(soil_nmds_df, aes(x = NMDS1, y = NMDS2, color = Type, fill = Type)) +
  geom_point(shape = 21, size = 1, stroke = 1) +
  geom_polygon(data = soil_ellipse_df, aes(x = NMDS1, y = NMDS2, fill = Type, group = Type, color = Type),
               alpha = 0.4, linewidth = 0.5) +
  # Species labels
  geom_text(data = soil_species_scores, aes(x = NMDS1, y = NMDS2, label = Species),
            color = "darkcyan", alpha = 0,
            inherit.aes = FALSE,
            size = 3) +
  scale_color_manual(name = "Forest Treatment",
                     values = colvec) +
  scale_fill_manual(name = "Forest Treatment",
                    values = colvec) +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
    strip.background = element_rect(fill = "grey", color = "black", size = 1),
    strip.text = element_text(color = "black", size = 12),
    panel.spacing = unit(1, "lines"),
    panel.background = element_rect(fill = "white", color = "black"),
    plot.margin = margin(10, 10, 10, 10)) +
  annotate("text", x = Inf, y = Inf, label = soil_anosim_label, hjust = 1.1, vjust = 1.5, size = 3.5, family = "mono") +
  labs(title = "Soil Invertebrate Community Separation")

soil_nmds_hull_graph + coord_fixed()

save_plot(soil_nmds_hull_graph, "soil_nmds_hull_graph")
## png 
##   2

Soil invertebrate NMDS resampling

library(vegan)

soil_invert_wide_nonzero <- soil_invert_wide[rowSums(soil_invert_wide) > 0, ]
soil_invert_env_nonzero <- soil_invert_env[rowSums(soil_invert_wide) > 0, ]


nonzero_rows <- rowSums(soil_invert_wide) > 0
soil_invert_wide_nonzero <- soil_invert_wide[nonzero_rows, ]
soil_invert_env_nonzero <- soil_invert_env[nonzero_rows, ]

anyNA(soil_invert_wide_nonzero)  # Should return FALSE
## [1] FALSE
library(vegan)

base_R <- anosim(
  wisconsin(sqrt(soil_invert_wide_nonzero)),
  soil_invert_env_nonzero$Type,
  distance = "bray",
  permutations = 999
)$statistic


# Use non-zero rows only
base_R <- anosim(wisconsin(sqrt(soil_invert_wide_nonzero)), soil_invert_env_nonzero$Type, distance = "bray", permutations = 999)$statistic


jackknife_results <- data.frame(Species = character(),
                                R_value = numeric(),
                                Delta_R = numeric(),
                                stringsAsFactors = FALSE)

species_names <- colnames(soil_invert_wide_nonzero)

for (sp in species_names) {
  # Remove the species (column)
  reduced_matrix <- soil_invert_wide_nonzero[, setdiff(species_names, sp)]
  
  # Remove rows that are now completely zero
  nonzero_rows <- rowSums(reduced_matrix) > 0
  reduced_matrix <- reduced_matrix[nonzero_rows, ]
  reduced_env <- soil_invert_env_nonzero[nonzero_rows, ]

  # Skip iteration if data becomes too small or empty
  if (nrow(reduced_matrix) < 2 || anyNA(reduced_matrix)) next
  
  # Run ANOSIM
  anosim_result <- anosim(wisconsin(sqrt(reduced_matrix)), reduced_env$Type, distance = "bray", permutations = 999)
  
  # Store results
  R_val <- anosim_result$statistic
  delta_R <- base_R - R_val
  jackknife_results <- rbind(jackknife_results, data.frame(Species = sp, R_value = R_val, Delta_R = delta_R))
}


# Sort
jackknife_results <- jackknife_results[order(-jackknife_results$Delta_R), ]

# Plot
soil_inv_jk_graph <- ggplot(jackknife_results, aes(x = reorder(Species, Delta_R), y = Delta_R)) +
  geom_col(fill = "#0072B2") +
  coord_flip() +
  labs(title = "Soil Invertebrate Order Impact on NMDS Clustering",
       x = "Order Removed",
       y = "Decrease in ANOSIM R") +
  theme_minimal()

soil_inv_jk_graph

save_plot(soil_inv_jk_graph, "soil_inv_jk_graph")
## png 
##   2

Ground Invertebrate Diversity:

Extracting model predicted ground invertebrate diversity:

predicted_pitfall_diversity <- predict(pitfall_diversity_model, type = "response", level = 1)

predicted_pitfall_div_with_ci <- intervals(pitfall_diversity_model, level = 0.95)

# Extract just the fixed effects into a data frame
pit_div_ci_df <- as.data.frame(predicted_pitfall_div_with_ci$fixed)


# Optionally, name the columns
colnames(pit_div_ci_df) <- c("CI_lower", "Estimate", "CI_upper")

# Add a column for the effect names (e.g., Intercept, Carbon, etc.)
pit_div_ci_df$Effect <- rownames(predicted_pitfall_div_with_ci$fixed)

# Reorder for readability
pit_div_ci_df <- pit_div_ci_df[, c("Effect", "Estimate", "CI_lower", "CI_upper")]

# View it
print(pit_div_ci_df)
##                                            Effect     Estimate    CI_lower
## (Intercept)                           (Intercept)  0.661053722  0.55147425
## TypeRegen                               TypeRegen -0.002566658 -0.12499009
## TypeUnforested (heath)     TypeUnforested (heath)  0.077685011 -0.05544383
## TypeUnforested (remnant) TypeUnforested (remnant) -0.004524457 -0.12391595
##                           CI_upper
## (Intercept)              0.7706332
## TypeRegen                0.1198568
## TypeUnforested (heath)   0.2108138
## TypeUnforested (remnant) 0.1148670
# Extract the intercept value
pit_intercept_estimate <- pit_div_ci_df$Estimate[pit_div_ci_df$Effect == "(Intercept)"]
pit_intercept_lower <- pit_div_ci_df$CI_lower[pit_div_ci_df$Effect == "(Intercept)"]
pit_intercept_upper <- pit_div_ci_df$CI_upper[pit_div_ci_df$Effect == "(Intercept)"]

# Now create a new column with actual predicted values
pit_div_ci_df$Actual_Estimate <- pit_div_ci_df$Estimate
pit_div_ci_df$Actual_CI_lower <- pit_div_ci_df$CI_lower
pit_div_ci_df$Actual_CI_upper <- pit_div_ci_df$CI_upper

# Add the intercept to all *non-intercept* rows
pit_non_intercept_rows <- pit_div_ci_df$Effect != "(Intercept)"

pit_div_ci_df$Actual_Estimate[pit_non_intercept_rows] <- 
  pit_div_ci_df$Estimate[pit_non_intercept_rows] + pit_intercept_estimate

pit_div_ci_df$Actual_CI_lower[pit_non_intercept_rows] <- 
  pit_div_ci_df$CI_lower[pit_non_intercept_rows] + pit_intercept_lower

pit_div_ci_df$Actual_CI_upper[pit_non_intercept_rows] <- 
  pit_div_ci_df$CI_upper[pit_non_intercept_rows] + pit_intercept_upper


# Using base R
colnames(pit_div_ci_df)[colnames(pit_div_ci_df) == "Effect"] <- "Type"

pit_div_ci_df$Type[pit_div_ci_df$Type == "(Intercept)"] <- "Mature"
pit_div_ci_df$Type[pit_div_ci_df$Type == "TypeRegen"] <- "Regen"
pit_div_ci_df$Type[pit_div_ci_df$Type == "TypeUnforested (heath)"] <- "Unforested (heath)"
pit_div_ci_df$Type[pit_div_ci_df$Type == "TypeUnforested (remnant)"] <- "Unforested (remnant)"

pit_div_ci_df
##                                          Type     Estimate    CI_lower
## (Intercept)                            Mature  0.661053722  0.55147425
## TypeRegen                               Regen -0.002566658 -0.12499009
## TypeUnforested (heath)     Unforested (heath)  0.077685011 -0.05544383
## TypeUnforested (remnant) Unforested (remnant) -0.004524457 -0.12391595
##                           CI_upper Actual_Estimate Actual_CI_lower
## (Intercept)              0.7706332       0.6610537       0.5514742
## TypeRegen                0.1198568       0.6584871       0.4264842
## TypeUnforested (heath)   0.2108138       0.7387387       0.4960304
## TypeUnforested (remnant) 0.1148670       0.6565293       0.4275583
##                          Actual_CI_upper
## (Intercept)                    0.7706332
## TypeRegen                      0.8904900
## TypeUnforested (heath)         0.9814470
## TypeUnforested (remnant)       0.8855002

Plotting ground invertebrate diversity

pitfall_diversity_graph <- ggplot(pitfall_div, aes(x = Type, y = Diversity_Index, fill = Type)) +

  # Violin plot for raw data distribution
  geom_violin(alpha = 0.3, color = NA, width = 0.9) +

  # Jittered individual data points
  geom_jitter(aes(color = Type), size = 2, shape = 21, width = 0.1, height = 0, alpha = 0.6) +

  # Model-estimated means
  geom_point(data = pit_div_ci_df,
             aes(x = Type, y = Actual_Estimate),
             shape = 21, fill = "black", color = "black", size = 3, inherit.aes = FALSE) +

  # Model-estimated confidence intervals
  geom_errorbar(data = pit_div_ci_df,
                aes(x = Type, ymin = Actual_CI_lower, ymax = Actual_CI_upper),
                width = 0.2, color = "black", size = 0.5, inherit.aes = FALSE) +

  # Clean, minimal theme
  theme_bw() +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1),
    strip.background = element_rect(fill = "grey", color = "black", size = 1),
    strip.text = element_text(color = "black", size = 12),
    panel.spacing = unit(1, "lines"),
    panel.background = element_rect(fill = "white", color = "black"),
    plot.margin = margin(10, 10, 10, 10),
    legend.position = "none"
  ) +
  labs(
    x = "Forest Treatment",
    y = "Diversity Index (1 - Simpson's D)",
    fill = "Type",
    title = "Ground Invertebrate Diversity"
  ) +
  scale_fill_manual(
    name = "Forest Treatment",
    values = c(
    "Mature" = "#D55E00",
    "Regen" = "#009E73",
    "Unforested (heath)" = "#56B4E9",
    "Unforested (remnant)" = "#CC79A7"
  )) +
  scale_color_manual(
    name = "Forest Treatment",
    values = c(
    "Mature" = "#D55E00",
    "Regen" = "#009E73",
    "Unforested (heath)" = "#56B4E9",
    "Unforested (remnant)" = "#CC79A7"
  ))

Pitfall Invertebrate Richness:

Extracting model predicted pitfall invertebrate richness:

# Step 1: Extract the intercept value
pit_intercept_value <- coef(summary(pitfall_richness_model))["(Intercept)", "Estimate"]

pit_wald_ci <- confint(pitfall_richness_model, level = 0.95)
## Computing profile confidence intervals ...
# Drop .sig01 row
pit_wald_ci <- pit_wald_ci[!rownames(pit_wald_ci) %in% ".sig01", ]

# Then build the dataframe
pit_results_df <- data.frame(
  Effect = rownames(coef(summary(pitfall_richness_model))),
  Estimate = coef(summary(pitfall_richness_model))[, "Estimate"],
  Std_Error = coef(summary(pitfall_richness_model))[, "Std. Error"],
  p_value = coef(summary(pitfall_richness_model))[, "Pr(>|z|)"],
  CI_Lower = pit_wald_ci[, 1],
  CI_Upper = pit_wald_ci[, 2]
)


# Step 4: Add real (absolute) estimates and CIs — only modify non-intercepts
pit_results_df$Real_Estimate <- pit_results_df$Estimate
pit_results_df$Real_CI_Lower <- pit_results_df$CI_Lower
pit_results_df$Real_CI_Upper <- pit_results_df$CI_Upper

pit_non_intercepts <- pit_results_df$Effect != "(Intercept)"

pit_results_df$Real_Estimate[pit_non_intercepts] <- pit_results_df$Estimate[pit_non_intercepts] + pit_intercept_value
pit_results_df$Real_CI_Lower[pit_non_intercepts] <- pit_results_df$CI_Lower[pit_non_intercepts] + pit_intercept_value
pit_results_df$Real_CI_Upper[pit_non_intercepts] <- pit_results_df$CI_Upper[pit_non_intercepts] + pit_intercept_value

# Step 5: Rename 'Effect' column to 'Type' and recode factor levels
colnames(pit_results_df)[colnames(pit_results_df) == "Effect"] <- "Type"

pit_results_df$Type[pit_results_df$Type == "(Intercept)"] <- "Mature"
pit_results_df$Type[pit_results_df$Type == "TypeRegen"] <- "Regen"
pit_results_df$Type[pit_results_df$Type == "TypeUnforested (heath)"] <- "Unforested (heath)"
pit_results_df$Type[pit_results_df$Type == "TypeUnforested (remnant)"] <- "Unforested (remnant)"

# Step 6: (Optional) Set row names to match the cleaned-up 'Type' values
rownames(pit_results_df) <- pit_results_df$Type



# Step 8: Exponentiate Estimates and Confidence Intervals
pit_results_df$Real_Estimate <- exp(pit_results_df$Real_Estimate)
pit_results_df$Real_CI_Lower <- exp(pit_results_df$Real_CI_Lower)
pit_results_df$Real_CI_Upper <- exp(pit_results_df$Real_CI_Upper)

print(pit_results_df)
##                                      Type   Estimate Std_Error      p_value
## Mature                             Mature  1.7717679 0.1278215 1.087372e-43
## Regen                               Regen -0.2433485 0.1347723 7.097659e-02
## Unforested (heath)     Unforested (heath) -0.1389162 0.1539754 3.669519e-01
## Unforested (remnant) Unforested (remnant) -0.1734879 0.1337063 1.944492e-01
##                        CI_Lower   CI_Upper Real_Estimate Real_CI_Lower
## Mature                1.4471914 2.08273351      5.881242      4.251158
## Regen                -0.5095904 0.02007054      4.610883      3.533106
## Unforested (heath)   -0.4462172 0.15916599      5.118450      3.764258
## Unforested (remnant) -0.4374303 0.08802427      4.944521      3.797480
##                      Real_CI_Upper
## Mature                    8.026379
## Regen                     6.000474
## Unforested (heath)        6.895947
## Unforested (remnant)      6.422402

Plotting ground invertebrate richness

pitfall_richness_graph <- ggplot(pit_results_df, aes(x = Type, y = Real_Estimate, fill = Type)) +
  geom_jitter(data = pitfall_rich, aes(x = Type, y = Richness, color = Type), 
              size = 2, shape = 21, width = 0.1, height = 0, alpha = 0.6) +  # Add individual data points
  geom_point(size = 3, shape = 21, color = "black", fill = "black", position = position_dodge(width = 0.8)) +  # Scatter plot with mean
  geom_errorbar(
    aes(ymin = Real_CI_Lower, ymax = Real_CI_Upper),  # Error bars based on confidence intervals
    width = 0.2,  # Error bar width
    color = "black",  # Error bar color
    size = 0.5  # Error bar size
  ) +
  geom_violin(data = pitfall_rich, aes(x = Type, y = Richness, fill = Type),
            alpha = 0.3, color = NA, width = 0.9) +
  theme_bw() +  # Minimal theme
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1),  # Rotate x-axis labels
    strip.background = element_rect(fill = "grey", color = "black", size = 1),  # Grey background for facet labels
    strip.text = element_text(color = "black", size = 12),
    panel.spacing = unit(1, "lines"),  # Add spacing between panels if needed
    panel.background = element_rect(fill = "white", color = "black"),  # Apply white background with black border to panels
    plot.margin = margin(10, 10, 10, 10)  # Text color inside the grey box
  ) +
  labs(
    x = "Forest Treatment",  # x-axis label
    y = "Order level Richness",  # y-axis label
    fill = "Type",
    title = "Ground Invertebrate Richness"  # Overall title
  ) +
  scale_fill_manual(
    name = "Forest Treatment",
    values = c(
    "Mature" = "#D55E00",
    "Regen" = "#009E73",
    "Unforested (heath)" = "#56B4E9",
    "Unforested (remnant)" = "#CC79A7"
  )) +
  scale_color_manual(
    name = "Forest Treatment",
    values = c(
    "Mature" = "#D55E00",
    "Regen" = "#009E73",
    "Unforested (heath)" = "#56B4E9",
    "Unforested (remnant)" = "#CC79A7"
  ))


pitfall_rich_div_graph <- pitfall_diversity_graph + pitfall_richness_graph

pitfall_rich_div_graph

save_plot(pitfall_rich_div_graph, "pitfall_rich_div_graph")
## png 
##   2

Pitfall Community Composition:

pitfall_dataset <- pitfall_dataset[c("Order", "Type", "Abundance", "Site.y", "Plot")]


pitfall_dataset <- aggregate(cbind(Abundance) ~ Order + Plot, data = pitfall_dataset, FUN = sum)

pitfall_dataset <- merge(pitfall_dataset, overview_dataset, by = "Plot")

pitfall_dataset <- pitfall_dataset[c("Order", "Type", "Abundance", "Site", "Plot")]

pitfall_wide <- as.data.frame(spread(pitfall_dataset, Order, Abundance))

After converting the data into wide format, I need to replace all the NA values with 0

pitfall_wide[is.na(pitfall_wide)] <- 0

pitfall_env <- pitfall_wide[c("Plot", "Type", "Site")]

pitfall_wide$Plot <- NULL
pitfall_wide$Site <- NULL
pitfall_wide$Type <- NULL
pitfall_wide$Family <- NULL
pitfall_wide$Genus <- NULL

pitfall_wide <- as.data.frame(lapply(pitfall_wide, as.integer))

str(pitfall_wide)
## 'data.frame':    74 obs. of  14 variables:
##  $ Acari           : int  0 0 1 1 0 0 1 0 1 3 ...
##  $ Araneae         : int  0 2 3 1 0 1 0 2 3 1 ...
##  $ Arionidae       : int  0 0 0 1 0 0 0 0 0 0 ...
##  $ Coleoptera      : int  0 0 2 6 13 1 3 1 7 6 ...
##  $ Diptera         : int  2 0 2 5 1 10 2 10 2 5 ...
##  $ Entomobryomorpha: int  2 0 3 1 2 7 1 6 5 14 ...
##  $ Hemiptera       : int  0 0 0 3 0 0 0 0 1 1 ...
##  $ Hymenoptera     : int  0 2 1 0 1 2 0 3 0 6 ...
##  $ Isopoda         : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Julida          : int  0 1 0 1 0 0 0 0 0 0 ...
##  $ Lepidoptera     : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Opisthopora     : int  0 0 1 1 0 0 1 0 0 0 ...
##  $ Podomorpha      : int  0 0 0 0 0 5 1 0 0 1 ...
##  $ Symphpleona     : int  0 0 0 0 0 0 0 0 0 0 ...
NMDS_pitfall <- metaMDS(pitfall_wide)
## Square root transformation
## Wisconsin double standardization
## Run 0 stress 0.2371 
## Run 1 stress 0.2403582 
## Run 2 stress 0.2379947 
## Run 3 stress 0.2325491 
## ... New best solution
## ... Procrustes: rmse 0.08307471  max resid 0.3866432 
## Run 4 stress 0.2327488 
## ... Procrustes: rmse 0.01792176  max resid 0.1000026 
## Run 5 stress 0.242124 
## Run 6 stress 0.2391028 
## Run 7 stress 0.2325653 
## ... Procrustes: rmse 0.01496153  max resid 0.1016296 
## Run 8 stress 0.2329172 
## ... Procrustes: rmse 0.0300439  max resid 0.1584433 
## Run 9 stress 0.2331662 
## Run 10 stress 0.2468438 
## Run 11 stress 0.233911 
## Run 12 stress 0.238065 
## Run 13 stress 0.2327124 
## ... Procrustes: rmse 0.04351231  max resid 0.2047893 
## Run 14 stress 0.2417877 
## Run 15 stress 0.2331679 
## Run 16 stress 0.2326182 
## ... Procrustes: rmse 0.01523301  max resid 0.1015965 
## Run 17 stress 0.2349177 
## Run 18 stress 0.2331387 
## Run 19 stress 0.2321276 
## ... New best solution
## ... Procrustes: rmse 0.02328379  max resid 0.103097 
## Run 20 stress 0.2624901 
## *** Best solution was not repeated -- monoMDS stopping criteria:
##      3: no. of iterations >= maxit
##     17: stress ratio > sratmax
NMDS_pitfall
## 
## Call:
## metaMDS(comm = pitfall_wide) 
## 
## global Multidimensional Scaling using monoMDS
## 
## Data:     wisconsin(sqrt(pitfall_wide)) 
## Distance: bray 
## 
## Dimensions: 2 
## Stress:     0.2321276 
## Stress type 1, weak ties
## Best solution was not repeated after 20 tries
## The best solution was from try 19 (random start)
## Scaling: centring, PC rotation, halfchange scaling 
## Species: expanded scores based on 'wisconsin(sqrt(pitfall_wide))'
stressplot(NMDS_pitfall)

pitfall_env
##    Plot                 Type         Site
## 1    63               Mature Glen Cluanie
## 2    66               Mature Glen Cluanie
## 3    68               Mature Glen Cluanie
## 4    69               Mature Glen Cluanie
## 5    74               Mature Glen Cluanie
## 6    81               Mature Glen Cluanie
## 7    82               Mature Glen Cluanie
## 8    36               Mature   Glen Loyne
## 9    39               Mature   Glen Loyne
## 10   40               Mature   Glen Loyne
## 11   41               Mature   Glen Loyne
## 12   43               Mature   Glen Loyne
## 13   48               Mature   Glen Loyne
## 14   50               Mature   Glen Loyne
## 15    3               Mature  Glen Mallie
## 16    4               Mature  Glen Mallie
## 17   10               Mature  Glen Mallie
## 18   12               Mature  Glen Mallie
## 19   14               Mature  Glen Mallie
## 20   24               Mature  Glen Mallie
## 21   25               Mature  Glen Mallie
## 22   62                Regen Glen Cluanie
## 23   64                Regen Glen Cluanie
## 24   65                Regen Glen Cluanie
## 25   67                Regen Glen Cluanie
## 26   71                Regen Glen Cluanie
## 27   73                Regen Glen Cluanie
## 28   80                Regen Glen Cluanie
## 29   35                Regen   Glen Loyne
## 30   37                Regen   Glen Loyne
## 31   38                Regen   Glen Loyne
## 32   42                Regen   Glen Loyne
## 33   44                Regen   Glen Loyne
## 34   45                Regen   Glen Loyne
## 35   47                Regen   Glen Loyne
## 36    6                Regen  Glen Mallie
## 37    7                Regen  Glen Mallie
## 38    8                Regen  Glen Mallie
## 39   11                Regen  Glen Mallie
## 40   13                Regen  Glen Mallie
## 41   23                Regen  Glen Mallie
## 42   26                Regen  Glen Mallie
## 43   57   Unforested (heath) Glen Cluanie
## 44   58   Unforested (heath) Glen Cluanie
## 45   59   Unforested (heath) Glen Cluanie
## 46   60   Unforested (heath) Glen Cluanie
## 47   77   Unforested (heath) Glen Cluanie
## 48   29   Unforested (heath)   Glen Loyne
## 49   30   Unforested (heath)   Glen Loyne
## 50   31   Unforested (heath)   Glen Loyne
## 51   53   Unforested (heath)   Glen Loyne
## 52   54   Unforested (heath)   Glen Loyne
## 53   16   Unforested (heath)  Glen Mallie
## 54   20   Unforested (heath)  Glen Mallie
## 55   61 Unforested (remnant) Glen Cluanie
## 56   70 Unforested (remnant) Glen Cluanie
## 57   72 Unforested (remnant) Glen Cluanie
## 58   75 Unforested (remnant) Glen Cluanie
## 59   76 Unforested (remnant) Glen Cluanie
## 60   78 Unforested (remnant) Glen Cluanie
## 61   79 Unforested (remnant) Glen Cluanie
## 62   32 Unforested (remnant)   Glen Loyne
## 63   33 Unforested (remnant)   Glen Loyne
## 64   34 Unforested (remnant)   Glen Loyne
## 65   46 Unforested (remnant)   Glen Loyne
## 66   49 Unforested (remnant)   Glen Loyne
## 67   51 Unforested (remnant)   Glen Loyne
## 68   52 Unforested (remnant)   Glen Loyne
## 69    1 Unforested (remnant)  Glen Mallie
## 70    2 Unforested (remnant)  Glen Mallie
## 71    5 Unforested (remnant)  Glen Mallie
## 72   15 Unforested (remnant)  Glen Mallie
## 73   22 Unforested (remnant)  Glen Mallie
## 74   27 Unforested (remnant)  Glen Mallie
pitfall_env$Type <- factor(pitfall_env$Type)



scl <- 4
colvec <- c("#D55E00", "#009E73", "#56B4E9", "#CC79A7")
plot(NMDS_pitfall, type = "n")
with(pitfall_env, points(NMDS_pitfall, display = "sites", col = colvec[Type],
                          scaling = scl, pch = 21, bg = colvec[Type]))
## Warning in plot.xy(xy.coords(x, y), type = type, ...): "scaling" is not a
## graphical parameter
with(pitfall_env, legend("topleft", legend = levels(Type), bty = "n",
                          col = colvec, pch = 21, pt.bg = colvec))

NMDS_pitfall.fit <- envfit(NMDS_pitfall ~ Type , data = pitfall_env, perm = 999)
NMDS_pitfall.fit
## 
## ***FACTORS:
## 
## Centroids:
##                            NMDS1   NMDS2
## TypeMature                0.0521 -0.1043
## TypeRegen                -0.0412 -0.0321
## TypeUnforested (heath)   -0.0783  0.3440
## TypeUnforested (remnant)  0.0355 -0.0632
## 
## Goodness of fit:
##          r2 Pr(>r)
## Type 0.0528  0.243
## Permutation: free
## Number of permutations: 999
plot(NMDS_pitfall.fit, dis = "site")
## Warning in text.default(x$factors$centroids[, choices, drop = FALSE], labs$f, :
## "dis" is not a graphical parameter
plot(NMDS_pitfall.fit)

### Ground invertebrate NMDS

# 1. Calculate ellipses using vegan::ordiellipse() but don't draw them
library(vegan)

pit_ordi_ellipses <- ordiellipse(NMDS_pitfall, pitfall_env$Type, kind = "se", conf = 0.95, draw = "none")

# 2. Extract ellipse coordinates into a data frame
pit_ellipse_df <- do.call(rbind, lapply(names(pit_ordi_ellipses), function(group) {
  pit_ellipse_points <- vegan:::veganCovEllipse(pit_ordi_ellipses[[group]]$cov, 
                                            pit_ordi_ellipses[[group]]$center, 
                                            pit_ordi_ellipses[[group]]$scale)
  data.frame(pit_ellipse_points, Type = group)
}))
colnames(pit_ellipse_df)[1:2] <- c("NMDS1", "NMDS2")


library(ggplot2)

# Your main NMDS points
pit_nmds_df <- data.frame(
  NMDS1 = NMDS_pitfall$points[, 1],
  NMDS2 = NMDS_pitfall$points[, 2],
  Type = pitfall_env$Type
)

pit_anosim_plant_com_r_tr <- anosim(wisconsin(sqrt(pitfall_wide)), pitfall_env$Type, distance = "bray", permutations = 9999)


pit_R_val <- round(pit_anosim_plant_com_r_tr$statistic, 3)
pit_p_val <- ifelse(pit_anosim_plant_com_r_tr$signif < 0.001, "< 0.001", paste0("= ", signif(pit_anosim_plant_com_r_tr$signif, 2)))

pit_anosim_label <- paste0("ANOSIM R = ", pit_R_val, "\n", "p ", pit_p_val)

# Site scores
pit_site_scores <- as.data.frame(NMDS_pitfall$points)

colnames(pit_site_scores)[1:2] <- c("NMDS1", "NMDS2")

pit_site_scores$Type <- pitfall_env$Type  # Add grouping

# Species scores (usually in $species, or use scores())
pit_species_scores <- as.data.frame(scores(NMDS_pitfall, display = "species"))
pit_species_scores$Species <- rownames(pit_species_scores)

pit_nmds_hull_graph <- ggplot(pit_nmds_df, aes(x = NMDS1, y = NMDS2, color = Type, fill = Type)) +
  geom_point(shape = 21, size = 1, stroke = 1) +
  geom_polygon(data = pit_ellipse_df, aes(x = NMDS1, y = NMDS2, fill = Type, group = Type, color = Type),
               alpha = 0.4, linewidth = 0.5) +
  geom_text(data = pit_species_scores, aes(x = NMDS1, y = NMDS2, label = Species),
            color = "darkcyan", alpha = 0,
            inherit.aes = FALSE,
            size = 2) +
  scale_color_manual(name = "Forest Treatment",
                     values = colvec) +
  scale_fill_manual(name = "Forest Treatment",
                    values = colvec) +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),  # Rotate x-axis labels
    strip.background = element_rect(fill = "grey", color = "black", size = 1),  # Grey background for facet labels
    strip.text = element_text(color = "black", size = 12),
    panel.spacing = unit(1, "lines"),  # Add spacing between panels if needed
    panel.background = element_rect(fill = "white", color = "black"),  # Apply white background with black border to panels
    plot.margin = margin(10, 10, 10, 10) ) +
  annotate("text", x = Inf, y = Inf, label = pit_anosim_label, hjust = 1.1, vjust = 1.5, size = 3.5, family = "mono") +
  labs(title = "Ground Invertebrate community separation")

pit_nmds_hull_graph + coord_fixed()

save_plot(pit_nmds_hull_graph, "pit_nmds_hull_graph")
## png 
##   2

Ground invertebrate resampling

# Start by calculating the base ANOSIM R-value (with all species included)
base_R_pitfall <- anosim(wisconsin(sqrt(pitfall_wide)), pitfall_env$Type, distance = "bray", permutations = 999)$statistic

# Initialize a data frame to store the jackknife results
jackknife_results_pitfall <- data.frame(Species = character(),
                                        R_value = numeric(),
                                        Delta_R = numeric(),
                                        stringsAsFactors = FALSE)

# Loop over all species (columns in your wide matrix)
species_names_pitfall <- colnames(pitfall_wide)

for (sp in species_names_pitfall) {
  # Remove the species (column)
  reduced_matrix_pitfall <- pitfall_wide[, setdiff(species_names_pitfall, sp)]
  
  # Remove rows that are now completely zero
  nonzero_rows_pitfall <- rowSums(reduced_matrix_pitfall) > 0
  reduced_matrix_pitfall <- reduced_matrix_pitfall[nonzero_rows_pitfall, ]
  reduced_env_pitfall <- pitfall_env[nonzero_rows_pitfall, ]
  
  # Skip iteration if data becomes too small or empty
  if (nrow(reduced_matrix_pitfall) < 2 || anyNA(reduced_matrix_pitfall)) next
  
  # Run ANOSIM on the reduced matrix
  anosim_result_pitfall <- anosim(wisconsin(sqrt(reduced_matrix_pitfall)), reduced_env_pitfall$Type, distance = "bray", permutations = 999)
  
  # Store R and change from base
  R_val_pitfall <- anosim_result_pitfall$statistic
  delta_R_pitfall <- base_R_pitfall - R_val_pitfall
  jackknife_results_pitfall <- rbind(jackknife_results_pitfall, data.frame(Species = sp, R_value = R_val_pitfall, Delta_R = delta_R_pitfall))
}

# Sort by greatest impact (largest Delta_R)
jackknife_results_pitfall <- jackknife_results_pitfall[order(-jackknife_results_pitfall$Delta_R), ]

# Plot Jackknife Results
pitfall_jk <- ggplot(jackknife_results_pitfall, aes(x = reorder(Species, Delta_R), y = Delta_R)) +
  geom_col(fill = "#0072B2") +
  coord_flip() +
  labs(title = "Ground invertebrate Order impact on NMDS Clustering",
       x = "Species Removed",
       y = "Decrease in ANOSIM R") +
  theme_minimal()


 pitfall_jk

 save_plot(pitfall_jk, "pitfall_jk")
## png 
##   2